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Final  Report 


RESEARCH  ON  SPECIALIZED 
COMPUTATIONAL  METHODS  FOR 
FLUID-STRUCTURE  INTERACTION 
SIMULATIONS  FOR  ADVANCED  SUBMARINE 

TECHNOLOGY 


1  Introductory  Comments 

This  document  summarizes  research  done  during  the  period  November  15,  1988-February 
2,  1992  on  a  project  aimed  at  the  development  of  new  advanced  computational  methods  for 
modeling  coupled  elastic  scattering  problems  in  linear  acoustics.  The  principal  feature  of  the 
research  was  to  develop  and  explore  the  use  of  advanced  /ip-adaptive  finite  element  methods 
for  this  class  of  problems,  including  the  possible  use  of  coupled  boundary  element  and  finite 
element  methods. 

Early  in  the  project,  the  research  team  elected  to  take  a  fresh  approach  to  these  types  of 
modeling  problems:  it  was  known  that  large  scale  problems  in  linear  acoustics  can  represent 
enormous  computational  efforts,  and  that  new  approaches  would  be  necessary  in  order  to 
cope  with  the  increasing  complexity  of  a'-oustical  analyses  needed  in  the  design  of  acoustical 
systems  associated  with  modern  submarines.  Previous  experience  and  successes  by  members 
of  the  research  team  in  developing  new  adaptive  hp  finite  element  methods  for  problems  in 
solid  mechanics  and  fluid  mechanics  naturally  led  to  the  plan  to  develop  new  hp  strategies 
for  structural  acoustics.  But  having  made  this  decision,  the  important  question  of  how  these 
physical  phenomena  are  to  be  modeled  remain  the  next  major  decision  in  determining  the 
direction  of  the  research. 

Certainly,  the  classical  approaches  to  this  problem  are  well  documented:  acoustical  fluid 
is  modeled  using  the  Helmholtz  equation,  the  structure  is  modeled  using  linear  structural 
mechanics,  and  mechanical  coupling  is  provided  at  the  interface.  But  there  are  many  different 
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techniques  that  can  be  used  to  model  these  types  of  problems.  Traditionally,  the  problem 
is  solved  in  the  frequency  domain  with  a  boundary  element  type  formulation  to  modeling 
acoustical  fluid.  There  are,  however,  a  number  of  other  approaches  worth  considering.  These 
include  the  use  of  finite  elements  to  model  the  fluid,  with  appropriate  boundary  conditions  in 
the  far  field.  Also,  which  particular  boundary  element  method  would  be  best  suited  for  this 
class  of  problems  was,  at  the  beginning  of  this  effort,  unknown.  Codes  exist  which  employ  one 
boundary  element  formulation  or  another,  but  all  of  these  were  known  to  possess  a  number  of 
serious  deficiencies,  and  many  of  these  precluded  the  use  of  high-order  adaptive  hp  methods. 
Then  there  was  also  a  classical  approach  in  which  linear  acoustics  was  modeled  by  integrating 
a  system  of  hyperbolic  conservation  laws.  This  is  by  no  means  a  common  approach  to  these 
problems,  but  it  carried  with  it  the  possibility  of  incorporating  additional  physics  into  the 
formulation  that  may  be  of  use  later.  Moreover,  it  was  not  known  at  the  onset  how  these 
more  direct  methods  compared  in  efficiency  and  accuracy  with  some  of  the  conventional 
methods.  Finally,  mathematical  and  convergence  properties  of  various  boundary  element 
methods  have  only  begun  to  be  established  and  a  posteriori  error  elements,  vital  for  the  use 
of  any  adaptive  scheme,  were  not  available  for  boundary  element  methods,  and  in  particular, 
were  not  available  forcoupled  finite  element  and  boundary  element  methods. 

After  a  detailed  study  of  the  literature  and  a  number  of  initial  calculations,  two  parallel 
approaches  were  attempted  in  which  the  principal  effort  was  directed  toward  a  coupled 
boundary  element /finite  element  method  formulation  and  a  smaller  pilot  effort  was  initiated 
on  the  use  of  /ip-finite  element  techniques  for  analyzing  transient  coupled  acoustics  problems. 
This  latter  approach  occupied  around  15  to  20  percent  of  the  total  effort  and  was  considered 
from  the  onset  as  a  completely  exploratory,  somewhat  high-risk  component  of  the  work. 

The  prinicpal  mathematical  model  that  was  selected  to  provide  the  basis  of  this  study 
was  that  of  the  classical  problem  of  elastic  scattering  in  linear  acoustics:  an  elastic  body 
occupying  a  bounded  domain  fi  is  submerged  into  an  infinite  acoustical  fluid  subjected  to 
an  incident  pressure  field.  One  wishes  to  determine: 

•  the  pressure  field  in  the  fluid,  in  particular  on  the  surface  of  the  scatterer,  and 

•  the  velocity  field  in  the  elastic  structure. 

The  time  variable  is  eliminated  using  Fourier  transforms  and  the  problem  reduces  the  solu¬ 
tion  to  the  Helmholtz  equation  in  the  exterior  domain,  coupled  with  the  equations  of  linear 
elasticity  in  ft  in  the  frequency.  The  principal  idea,  therefore,  was  similar  to  that  in  the 
NASHUA  code  [2].  The  boundary  integral  approach  is  used  to  replace  the  Helmholtz  equa¬ 
tion  and  the  problem  reduces  to  the  solution  of  linear  elasticity/lincar  structural  mechanics 
equations  inside  H  coupled  with  a  boundary  integral  formulation  on  the  boundary  F  =  dQ. 


In  both  the  BEM/FEM  study  and  in  the  transient  acoustic  study,  new  data  structures 
were  developed  that  would  allow  the  implementation  of  a  new  class  of  adaptive  hp  finite 
element  methods  in  which  both  the  local  mesh  size  h  and  the  order  p  of  the  approximation 
were  treated  as  free  parameters.  In  addition,  it  is  necessary  to  develop  mathematically  rig¬ 
orous  a  posteriori  error  estimates  to  drive  the  /ip-adaptive  process.  The  ultimate  purpose  of 
developing  such  high  order  adaptive  schemes  was  an  effort  to  produce  exponential  rates  of 
convergence,  so  that,  it  was  hoped,  highly  accurate  results  could  be  obtained  with  relatively 
few  degrees  of  freedom  compared  with  conventional  methods.  This  super-algebraic  conver¬ 
gence  property,  which  was  proved  earlier  to  hold  for  linear  elliptic  problems,  if  also  applicable 
to  the  subject  classes  of  problems  would  allow  one  of  the  principal  missions  of  the  project  to 
be  fulfilled:  the  treatment  of  very  large-scale  modeling  problems  with  relatively  few  degrees 
of  freedom.  In  summary,  the  principal  goals  of  the  project  can  be  stated  as  follows: 

•  The  derivation  of  a  uniformly  stable  boundary  element  formulations  in  the  equations 
valid  for  the  whole  range  of  wave  numbers  k.  The  formulation  should  be  compati¬ 
ble  with  the  the  regularity  assumption  for  finite  element  approximations  of  the 
structure. 

•  The  development  of  (7°  finite  element  approximations  for  one-,  two-,  and  three- 
dimensional  geometries.  This  part  of  the  project  included  such  fundamental  issues 
as  the  development  of  new  hp  finite  element  data  structures  and  constrained  approxi¬ 
mations  for  coupled  boundary  element /finite  element  problems. 

•  The  development  of  a  coupled  finite  element /boundary  element  code  and  the  veri¬ 
fication  of  expected  rates  of  convergence.  This  included  first  the  development  of  a 
two-dimensional  code  and  led  to  preliminary  work  on  a  three-dimensional  code  that 
is  currently  under  study  and  further  development. 

•  The  derivation,  implementation,  and  verification  of  an  a  posteriori  error  estimate  for 
boundary  element  methods.  These  estimates  perhaps  represent  a  significant  variance 
from  those  that  have  been  reported  for  the  simplest  one-dimensional  cases. 

•  The  development  of  adaptive  strategies.  Having  obtained  a  local  error  estimate,  these 
strategies  would  provide  for  reducing  mesh  size  or  increasing  spectral  order  so  as  to 
achieve  a  target  error  with  a  number  of  unknowns. 

•  The  solution  of  selected  problems.  These  are  designed  to  demonstrate  the  effective¬ 
ness  of  the  methodology  on  classical  benchmark  problems  and,  ultimately,  on  realistic 
scattering  problems  of  interest  in  Naval  research. 
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•  The  verification  of  the  applicability  of  the  hp  approximation  to  structural  mechanics, 
including  bulky  structures  such  as  bulkheads,  but  also  possibly  shells,  thin  beams  and 
plates. 

•  The  development  of  a  mathematical  foundation  for  transient  acoustics,  including  study¬ 
ing  mathematical  properties  of  linear  hyperbolic  systems  associated  with  transient 
acoustics,  the  development  of  new  transient  fep-adaptive  schemes,  the  development  of 
high  order  transient  solvers  and  associated  a  posteriori  error  estimates,  and  associated 
adaptive  strategies. 

2  State-of-the-Art  at  the  Start 

It  is  accurate  to  say  that  at  the  beginning  of  this  project  virtually  no  papers  existed  in  the 
literature  on  large-scale  transient  acoustics  calculations.  The  theory  itself,  while  existing 
within  the  framework  of  abstract  Cauchy  problems  in  functional  analysis,  is  not  generally 
known  in  the  acoustics  literature  nor  were  the  investigators  able  to  find  large  scale  projects 
of  any  kind  that  used  finite  element  methods  for  these  types  of  problems.  Also,  it  was  not 
known  whether  or  not  these  approaches  would  be  fruitful,  or,  if  they  worked,  if  they  would  be 
competitive  with  more  traditional  schemes  or  if  they  would  produce  new  insight  into  acous¬ 
tical  scattering  phenomena.  These  questions,  to  an  extent,  remain  open  but  considerable 
progress  was  made  in  developing  these  new  methodologies. 

For  the  boundary  integral  formulation  and  the  eissociated  boundary  element  methods, 
the  following  main  schools  of  research  were  identified: 

•  The  German  school,  represented  mostly  by  mathematicians  (Wendland,  Stefan,  Cos- 
tubal,  and  others). 

•  The  American  school,  represented  by  full-spectrum  from  theoretical  mathematics  to 
practical  engineering  calculations  (Hsiao  Kleinman,  Prelite,  Cruze,  Bannerjee,  Rizzo, 
and  others). 

•  The  British  school  of  acoustics  (Burton,  Miller,  Jones,  Unsell  and  others)  and  BEM 
(Brebbia,  etc.). 

•  The  French  school  of  FEMs  (Nedelec,  Hamoli,  and  others). 

•  Other  works  (such  as  the  work  of  Johnston  in  Sweden,  among  others). 

A  variety  of  different  boundary  integral  formhlations  w’ere  known  at  the  start  of  this  project 
and  some  mathematical  basis  for  them  haye  been  established,  but  the  bulk  of  the  litera¬ 
ture  which  focused  on  acoustics  and  finite  element  approximations  emplo3’ed  the  classical 
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collocation  approach.  There  are  papers  in  the  literature  on  adaptive  boundary  integral  for¬ 
mulations,  and  among  these  we  mention  the  work  of  Rencis,  Rank,  and  Stephen  but  very 
few  numerical  solutions  to  any  caises  other  than  simple  test  problems  were  found. 


3  Technical  Approach  at  the  Start 

First,  the  initial  approaches  to  the  BEM/FEM  formulation  are  discussed  as  these  comprise 
the  major  portion  of  the  initial  work.  Due  to  the  complexity  of  the  problem  and  many 
uncertain  open  questions  it  was  decided  that  a  simplified  class  of  models  would  be  first  for¬ 
mulated  which  minimized  the  complexity  of  coping  with  the  behavior  of  the  elastic  structure. 
This  idea  consists  of  replacing  the  effects  of  the  structure,  particular  the  non-local  Green’s 
operator,  with  a  local  spring-like  approximation,  reducing  the  problem  to  the  solution  of  a 
Helmholtz  equation  in  the  exterior  domain  with  a  Robin  (impedance)  boundary  condition 
on  the  fluid/structure  interface.  Three  features  of  this  approach  should  be  noted: 

1.  A  mathematical  structure  of  this  simplified  model  resembles  the  general  problem;  the 
same  function  spaces  are  used,  natural  eigenfrequencies  exist,  and  many  of  the  com¬ 
putational  issues  are  the  same. 

2.  Wherever  a  boundary  integral  formulation  is  used,  the  resulting  formulation  involves 
exactly  the  same  integrals  as  for  the  general  coupled  problem. 

3.  The  boundary  element  code  being  developed  for  the  simplified  model  could,  with  some 
modifications,  could  be  reconfigured  to  be  used  as  a  major  part  of  the  final  code  for 
the  general  coupled  problem. 

Emphasis  was  first  placed  on  the  two-dimensional  problem.  While  most  of  the  mathemat¬ 
ical  details  remain  the  same  in  three  dimensions,  associated  technical  details,  especially  those 
involving  /ip-adaptivity  are  much  simpler  and  easier  to  implement  for  the  two-dimensional 
case.  Moreover,  it  was  unclear  in  the  beginning  of  the  project  precisely  what  type  of  bound¬ 
ary  formulation  would  be  most  appropriate.  It  was  known  that  the  traditional  Helmholtz 
integral  equation  or  the  hypersingular  integral  formulation,  or  the  more  popular  Burton- 
Miller  formulation  all  had  mathematical  and  technical  deficiencies  that  had  to  be  coped 
with.  Initial  calculations  were  nevertheless  done  with  the  Helmholtz  formulation  and  with  a 
form  of  the  hypersingular  integral  formulation. 

Secondly,  for  the  work  on  the  transient  acoustics,  a  number  of  essentially  one-dimensional 
cases  were  initially  explored.  A  few  of  the  issues  were  more  basic.  They  involved  a  study 
of  abstract  Cauchy  problems  for  linear  acoustics,  ignoring  initially  the  effects  of  the  elas¬ 
tic  structure.  It  also  involved  the  development  of  high  order  temporal  integration  schemes 
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to  advance  the  solution  in  time  with  an  accuracy  commiserate  with  the  order  of  the  spa¬ 
tial  approximation.  It  is  natural,  that  some  of  the  more  classical  high  order  schemes  were 
first  explored  including,  for  example,  higher  order  multistep  methods  such  as  the  Adams- 
Bashford  scheme,  various  types  of  Runge-Kutta  methods,  including  explicit  Runge-Kutta, 
implicit  Runge-Kutta,  and  singly  implicit  Runge-Kutta  methods,  and  some  initial  studies  of 
Taylor-Galerkin  schemes  which  heretofore  had  not  been  applied  to  these  types  of  acoustics 
formulations. 


4  Progress,  Success  and  Failures 

The  first  stage  of  this  project  was  completed  in  the  summer  of  1990  and  documented  in  [3]. 

At  this  point  the  following  decisions  were  made  and  results  were  obtained: 

•  The  existence  of  fictitious  frequencies  in  both  Helmholtz  formulation  and  the  standard 
hypersingular  formulation  were  found  to  be  major  deficiencies  in  these  techniques,  even 
though  they  provide  the  basis  of  a  number  of  commercial  and  laboratory  codes  that 
are  in  use  today.  The  Burton-Miller  formulation,  while  overcoming  this  particular 
deficiency,  nevertheless  had  also  a  number  of  deficiencies  which  are  not  adequately 
addressed  in  the  literature.  These  included  complexities  with  evaluating  hypersingular 
integrals,  and  with  the  overall  robustness  of  the  scheme  when  applied  to  reasonable 
sample  problems. 

•  The  Burton-Miller  formulation  was  nevertheless  finally  selected  as  a  basis  for  the  FEM 
code.  However,  a  Galerkin  approach  rather  than  a  collocation  approach  was  selected. 

•  An  T^-residual  technique  was  developed  for  a  posteriori evrov  estimation.  It  was  derived 
mathematically,  verified  experimentally,  and  proved  to  be  asymptotically  exact  for  a 
significant  class  of  BEM  formulations. 

•  A  two-dimensional  experimental  hp  boundary  element  code  was  developed  for  the 
model  problem  discussed  in  the  previous  section.  This  code  was  cornpleted,  and  used 
to  solve  several  trial  problems  including  the  infinite  cylinder  problem  in  two  dimensions. 

•  Simple  A-adaptive  strategies  and  p-adaptive  strategies  based  on  the  residual  error  es¬ 
timate  and  on  the  equidistribution  principle  were  developed.  These  were  applied  to 
model  problems  and  became  functional  in  mid-1990. 

•  A  singly-implicit  Runge-Kutta  scheme  was  devised  to  be  used  in  conjunction  with  an 
/ip  method  and  a  finite  difference  method  for  temporal  approximations  in  the  transient 
acoustics  studies.  Preliminary  a  posteriori  error  estimates  were  derived,  these  being 
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based  on  residual  methods  for  the  spatial  domain  and  a  rather  standard  Runge-Kutta 
scheme  for  controlling  the  time  step. 

The  principal  difficulties  encountered  at  this  point  in  the  project  were  connected  with 
the  correct  implementation  of  the  hypersingular  part  of  the  Burton-Miller  formulation.  In 
particular,  initial  studies  showed  that,  contrary  to  previous  expectations,  nonexponential 
convergence  rates  were  obtained.  Nevertheless,  the  Burton-Miller  formulation  for  the  bound¬ 
ary  integral  methods  and  the  h-p  adaptive  schemes  for  the  structural  approximation  were 
regarded  as  viable  approaches  it  was  decided  to  go  forward  with  the  implementation  of  these 
methodologies  for  a  broader  class  of  problems.  Work  was  then  directed  in  three  directions: 

1.  Development  of  a  two-dimensional  code  for  the  full  coupled  problem, 

2.  development  of  a  three-dimensional  code  for  the  model  problem  with  simplified  Robin 
boundary  conditions,  and 

3.  the  exploration  of  new  temporal  schemes,  in  addition  to  Runge-Kutta  methods,  for 
handling  high-order  transient  acoustic  calculations. 

All  of  these  tasks  were  completed  by  the  end  of  1991  and  reported  in  [4]  and  [5].  At  this 
point  in  the  project,  a  number  of  the  major  conceptual,  mathematical  and  modeling  issues 
were  thought  to  be  understood  and  the  basis  for  modeling  a  considerably  more  general  class 
of  problems  was  established.  In  particular: 


•  The  issue  of  the  correct  implementation  of  the  hypersingular  operators  was  success¬ 
fully  solved  in  [5]  and  in  more  detail  later  in  [6].  This  involved,  as  noted  earlier,  a 
full  Galerkin  approximation  of  the  hypersingular  integral  formulation;  we  believe  this 
was  the  first  such  formulation  attempted  in  the  literature.  Fortunately,  subsequent 
calculations  showed  that  this  led  to  optimal  rates  of  convergence  so  that,  once  again, 
the  possibility  of  obtaining  exponentially  convergent  techniques  for  general  classes  of 
acoustics  problems  emerged. 

•  A  three-dimensional  h-p  boundary  element  code  for  the  model  problem  was  completed 
and  verified  by  solving  the  problem  of  elastic  scattering  on  a  sphere  and  on  a  finite 
cylinder. 

•  A  technique  for  obtaining  a  posteriori  error  estimates  was  coded,  implemented,  and 
tested  on  a  number  of  sample  problems. 

•  A  two-dimensional  //p-adaptive  coupled  boundary  element/finite  element  code  for  cou¬ 
pled  scattering  problems  was  completed.  This  represented  the  completion  of  a  major 
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research  tool  in  which  a  variety  of  computational  issues  could  be  studied.  Initially  this 
code  was  used  to  study  scattering  of  plane  waves  on  infinite  elastic  cylinders. 

•  The  feasibility  of  appllying  hp  approximations  for  structural  problems  was  tested  on 
a  cylindrical  elastic  shell.  The  effect  of  shell  thickness  on  the  behavior  of  both  the 
structure  and  on  scattering  was  explored.  It  was  discovered  that  for  simple  geometries 
such  as  this,  a  two-dimensional  finite  element  code  based  on  /ip-adaptive  methods  was 
perfectly  adequate  and  could  approximate  the  shell  behavior  for  thicknesses  ranging 
from  very  thin  shells  to  very  thick  shells.  On  the  other  hand,  some  deterioration  in  the 
quality  of  solutions  for  problems  with  corners  and  discontinuities  was  observed.  It  was 
clear  from  these  examples  that  in  future  calculations  the  ability  to  use  /i-adaptivity  in 
the  neighborhood  of  scatterers  such  as  corners  represents  a  vital  and  useful  part  of  the 
overall  strategy. 

•  A  new  class  of  high-order  Taylor-Galerkin  schemes  were  developed  and  proved  to 
be  unconditionally  stable  for  the  transient  acoustics  problem.  This  was  regarded  as  a 
significant  advance  in  existing  technology.  These  schemes  proved  to  be  very  robust  and 
accurate  and  outperformed  classical  Runge-Kutta  methods  on  a  significant  number 
of  test  problems.  Some  preliminary  adaptive  strategies  for  transient  acoustics  were 
studied  and  a  research  code  was  developed  that  produced  a  number  of  interesting 
results.  Still,  much  work  remains  to  be  done  in  this  area.  Current  results  appear  to 
be  sensitive  to  the  particular  adaptive  strategy  used,  and  further  work  is  needed  in 
developing  robust  and  precise  error  estimates  for  these  classes  of  problems. 

•  Also,  progress  was  made  in  developing  adaptive  techniques  for  the  transient  response 
of  elastic  structures  coupled  to  the  transient  acoustic  package.  These  results  also  look 
promising,  although  they  are  still  very  preliminary  in  nature. 

It  can  be  said  that  at  the  conclusion  of  this  project  a  number  of  major  new  results  have 
been  established  and  a  good  starting  point  for  the  continuation  of  the  work  into  three- 
dimensional  simulations  is  in  place.  Further  work  will  be  augmented  by  studies  on  parallel 
computation,  for  it  is  believed  that  with  the  use  of  parallel  algorithms,  additional  efficiencies 
can  be  obtained  which  enhance  further  those  expected  to  be  produced  by  the  /ip-adaptive 
algorithms. 

5  Final  Technical  Approach  and  Final  Goals 

As  a  summary  of  some  of  the  ideas  outlined  earlier,  the  following  assumptions  have  been 
established  for  current  and  future  work.  The  weak  or  variational  form  of  the  Btirton-Miller 


formulation,  coupled  with  the  residual  v?tiational  formulation  for  the  elasticity  equations, 
will  be  used  as  a  basis  for  further  work. 

•  General  h-p  approximations  of  both  boundary  and  domain  variational  formulations 
will  be  used. 

•  A  posteriori  estimates  will  be  derived  based  on  error  residual  methods  for  the  structuic 
and  estimates  for  the  boundary  element  method. 

•  Continued  work  will  be  done  on  /ip-adaptive  data  structures  and  adaptive  strategies 
which  could  lead  to  exponential  rates  of  convergence. 

•  Work  will  continue  on  modeling  coupling  the  transient  behavior  of  the  acoustical  fluid 
to  that  of  an  elastic  structure  and  with  the  high-order  multistage  Taylor- Galerkin 
methods  will  be  the  driving  algorithm  for  the  temporal  approximation. 

•  Studies  of  parallel  algorithms  will  be  undertaken. 

6  Problems  Perceived  and  Future  Directions 

The  work  on  the  final  three-dimensional  code  for  the  coupling  problem  initiated  in  November 
of  1991  remains  to  be  a  major  task  in  the  current  research.  The  following  tasks  should  be 
listed: 

•  Development  of  a  three-dimensional  FE  code  for  the  elasticity  equations  coupled  to 
the  boundary  element  method. 

•  A  study  of  techniques  for  using  /ip-FEM’s  for  treating  thin  flexible  structures  such  as 
beams,  plates,  and  shells. 

•  A  geometrical  CAD-like  code  is  needed  which  is  capable  of  modeling  solid  objects  and 
general  surfaces.  This  requires  the  development  of  a  new  preprocessor  including  an 
h-p  mesh  generator  for  nontrivial  geometries. 

•  Development  of  an  a  posteriori  estimate  for  the  three-dimensional  coupled  problem 
together  with  some  new  /ip-adaptive  strategies. 

•  Extension  of  the  transient  acoustic  work  to  fully  coupled  problems  and  then  to  three 
dimensions. 

•  Exploration  of  domain  decomposition  methods  and  coarse-grain  parallel  algorithms 
lor  the  solution  of  some  large  scale  problems  on  multiprocessor  architectures. 
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As  for  difficulties  yet  to  be  resolved,  on  the  theoretical  side  the  integration  of  almost 
singular  functions  remains  a  critical  issue.  The  development  of  inexpensive  but  precise 
integration  procedures  for  almost  singular  integrands  may  determine  the  overall  success  of 
the  proposed  methodology.  Currently  rather  simple  adaptive  procedures  are  used  and  these 
may  prove  to  be  quite  expensive  unless  more  simplified  quadratures  can  be  obtained.  On  the 
other  hand,  simplified  quadrature  rules  may  defeat  all  the  gains  obtained  using  adaptivity. 

Research  on  a  theory  of  a  posteriori  error  estimation  must  be  continued.  This  includes 
the  study  of  more  practical  inexpensive  implementations  of  residual  estimates  and  a  search 
for  new  estimates  that  would  be  more  sensitive  to  the  errors  in  derivatives  in  the  solution. 
Most  importantly,  techniques  for  handling  geometric  singularities  such  as  lines,  corners  and 
material  interfaces  which  are  known  to  be  very  important  in  the  scattering  of  acoustical 
waves  must  be  studied.  The  major  issues  have  to  do  with  the  geometry  of  these  classes  of 
problems  and  how  well  that  geometry  is  approximated.  Our  feeling  is  that  the  geometry 
questions  are  open,  and  that  a  considerable  amount  of  additional  work  is  needed  before 
they  are  adequately  understood  or  correctly  factored  into  any  reasonable  acoustical  analysis 
capability. 
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Appendix  A 


* 


Variational  Formulations  and  /ip- Boundary 
Element  Approximations  for  Hyper  singular  Integral 
Equations  for  Helmholtz  Exterior  Boundary-Value 
Problems  in  Two  Dimensions 


Abstract 

In  this  paper,  a  weak  hypersingular  formulation  of  the  Helmholtz  exterior  bound¬ 
ary-value  problem  in  two  dimensions  is  presented.  The  weak  formulation  derived 
here  is  implemented  into  an  h-p-adaptive  boundary  element  approximation  and  the 
elementwise  L^-residual  is  used  as  the  error  estimation  in  the  adaptive  scheme.  A 
series  of  numerical  experiments  on  the  convergence  of  the  solution  and  the  properties 
of  T^-residual  is  given.  The  work  is  an  extension  of  methods  and  results  contained  in 
[2]  and  [3]. 

1  Introduction 

The  aim  of  this  paper  is  to  formulate  weak  hypersingular  boundary  integral  equations  for 
the  exterior  boundary-value  problems  for  the  Helmholtz  equations  in  two  dimensions.  The 
presented  work  is  an  extension  of  methods  and  results  of  earlier  work  [2].  Additional  details 
concerning  formulation  of  the  boundary-value  problem,  its  physical  meaning  and  related 
mathematical  concepts  can  be  found  in  [2]  and  [3]. 

The  priority  in  the  development  of  variational  formulations  of  hypersingular  integral  equa¬ 
tions  for  the  Helmholtz  exterior  boundary-value  problem  belongs  to  M.  A.  Hamdi  [5].  The 

*On  leave  from  the  Section  of  Applied  Mathematics  of  Technical  University  of  Cracow,  31-55  Krakow, 
ul.  Warszawska  24,  Poland. 
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classical  Helmholtz  exterior  hoiindary  value  problem  is  the  following:  For  a  given  bounded 
domain  Q  with  a  smooth  boundary  1'  and  completion  il"  (Fig.  1),  (ind  a  function  —*  C 
which  satisfies 

•  the  Helmholtz  differential  equation 

—  Au  —  k^u  =  0  in  H', 

•  the  Sommerfeld  boundary  condition  at  infinity 

=  o(r-5)  (1.2) 


(1.1) 


•  a  boundary  condition  on  F, 


where  k  is  the  wave  number. 


In  a  mathematical  description  of  acoustic  scattering,  a  scattered  pressure  p‘  is  cissumed 
to  fulfill  equations  (1.1)  and  (1.2).  The  incident  pressure  fulfills  equation  (1.1)  in  the 

whole  space  and  the  total  pressure 

p  =  f 

(1.3) 

fulfills  one  of  the  typical  boundary  conditions  on  F: 

p  =  0  or 

(1.4) 

9p 

—  =  0  or 

UTl^ 

(1.5) 

dp 

IhT,  = 

(1.6) 

In  ( l.fi),  £  is  a  positive  constant.  (1  f  =  0,  w<'  obtain  the  for 

inulation  for  rigid  ''Cattering;  for 

£  >  0,  a  formulation  of  a  compliant  boundary  with  spring  stiffness  is  obtained 


2  Classical  Integral  Formulations 

It  is  well  known  that  the  funda  nental  solution  of  (1.1)  in  two  dimensions  has  a  form  (cf.  [1]) 

=  '-Ho(kr)  (2.1) 
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where  =  -1,  r  =  ®  -  y,  r  =  |r|,  and  Hq  is  the  Hankel  function  of  the  first  kind  and  zeroth 
order.  The  function  <f>  may  also  be  written  as 

2/)  =  ^  +  ^o(r)  (2.2) 

where  ^o(t')  €  C°°(0,oo)  and  <^o  with  its  derivative  are  bounded  in  the  neighborhood  of  0. 
We  assume  next  that 

uec^  (O')  n  (n‘)  (2.3) 

where  C*’"  is  the  space  of  functions  the  first  derivatives  of  which  fulfill  a  Holder  condition 
with  exponent  q. 

From  the  Green  formula,  taking  into  account  (1.1)  and  (1.2)  we  may  obtain  the  Helmholtz 
Integral  Formulation  on  the  boundary  F  ([3]): 

^u{x)  =  ^\u{y)-^{x,y)  - -^{y)(f>{x,y)  ds(y) ,  x  £  T  (2.4) 

L  2/ 

where  the  integral  is  understood  as  a  Cauchy  Principal  Value  Integral,  however  it  has  been 
proven  [3]  that  it  is  in  fact  a  Lebesgue  integral. 

For  the  pressure,  the  Helmholtz  Integral  Formulation  becomes 

^p{x)  =  J^  ^{x,y)p{y)  -  <f>{x,y)-^{y)  ds{y)  +  p'’'‘ {x)  (2.5) 

•  y  y  • 

because  multiplied  by  <{>  satisfies  the  second  Green  identity  inside  Q: 

=  o  (2.6) 

y  y 

Using  a  similar  technique,  the  hypersingular  integral  formulation  may  be  obtained  [3] 


which,  for  the  pressure,  becomes 

(2.8) 

The  first  integrals  in  both  equations  should  be  considered  as  Hadainard  Finite  Part  Integrals; 


Equations  (2.4)  and  (2.7)  may  be  obtained  as  limit  cases  of  corresponding  equations 
taken  for 


uM  =  I 


y 


_  dni 


ds{y) 


(2.10) 


and 


J^u{y)^^^(x„,y)ds{v)- (2.11) 

using  appropriate  limit  theorems  for  potentials.  Here  denotes  an  outward  normal  vector 
to  the  boundary  F  at  x  =  Jii^Xm  G  F.  For  the  total  pressure  p,  equations  (2.10)  and  (2.11) 
become 


d^(t> 


P(*m) 


-H 


3<f>  d 


y 


dn 


y 


d5(y)  +p’"‘=(x„) 


(2.12) 


dp 

dnx 


(*m) 


=x[ 


/  N  /  X  dp  d(t>  ,  ^ 

p(y)a„  •  -\^m,y)  -  ^(y)a:^(®m,y) 


dnx'dny 


drii 


dp' 


ds{y)^-^{xJi  (2.13) 


The  integrals  (2.10)-(2.13)  are  usual  Lebesgue  integrals  of  bounded  functions. 


3  Propositions  and  Lemmas 

Lemma  1.  If  ^(x,  y)  is  the  fundamental  solution  of  (1.1)  defined  by  (2.1),  then  the  following 
equality  holds: 

k'^  I  <f){x,y)dy+  [  ■^{x,y)ds{y)  +  l-  =  0,  VxGF  (3.1) 

Jii  Jr  o^y  J- 

Proof:  Let  us  define  a  ball 

BM  =  {y  ^S?-\\x-y\\<e]  (3.2) 

with  its  boundary  (circle) 

•?e(®)  =  {y  -  y||  =  e}  (3.3) 

(see  Fig.  2).  Applying  the  first  Green  identity  to  u{x)  =  1  and  (j>  we  obtain 

-I  ^y<l>{x,y)dy+  f  ■^{x,y)ds{y)+  I  ^^(x,y)ds(y)  =  0  (3.4) 

jn'\e,  Jr\B,  ony  Js.nn'  dny 

<t>  fulfills  (1.1),  then 

“  /  ^y<Pi^^y)dy  =  I  <i>{x,y)dy  -*  f  (f>{x,y)dy  (3.5) 
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On  the  circle  St  the  relation 


^(x,y)  =  -<A'(r) 

is  valid  (Fig.  2).  Then  for  the  smooth  boundary 

=  [5^ + ^'(^)]  X.„„.  ^  0  =  5 

The  second  integral  tends  to 

which,  with  (3.6),  gives  (3.1).  I 

Proposition  1.  For  a  piecewise  smooth  boundary  F,  any  p  G  C*(r)  and  any  ip(x,y)  = 
<p(r)  6  <^^(0, 00),  we  have 


where  ry  =  (y)  is  the  tangent  vector  at  y. 

Proof:  It  is  well  known  that  for  each  closed,  piecewise  smooth  curve  F  C  and  for  any 
function  u  €  C'^(r) 

^  Vu(y)  •  T{y)ds{y)  =  0  (3.9) 

If  we  substitute  now 

«(y)  =  ^(®.y)p(y)  (3.io) 

we  obtain  (3.8).  I 

Proposition  2.  For  any  function  (f  €  C'*(fi)  the  following  equality  holds: 

d  ( d^p  \ 

Proof:  We  use  formulas  (see  Fig.  3) 


r\ 

T2 


rij 

-n, 


(3.12) 


+  r.T,Ayv  =  _„.(x)„.(y)  +  _„,(xK(y) 

o2  ^2  ii2  cs2 

=  Ti(®)Ti(y)|-^  +  Ti(a!)r2(y)^^-^  +  r2(g)ri(y)-^  ^  +  r2(x)r2(y)-^ 


dyidya 


dy2dyi 


dyl 


=  T.(y) 


5v?  ,  J  d  f  dip  \  d  (  dip  \ 

^  dyi  ^  V^Ta;  j  5ry  (st®  j 


Let  us  fix  any  point  x  G  T  and  introduce  there  a  local  tangential-normal  coordinate 
system  (^1,^2)  (Fig.  3).  In  this  system,  the  boundary  F  will  be  given  by  a  function  / 


(616)  €  r  ^2  = /(^i) ,  V^i€  (62,^2)  <  0  <  62 


(3.13) 


Lemma  2.  Let  €  FI®,  Xm  — >  x  €  F.  Let  the  function  /  given  by  (3.13)  belong  to 
C*’"(^i,^2)  ‘‘•nd  p  be  bounded  on  F.  Then 


J^H^m,y)p(y]ds(y}  — *  J^<^ix,y]p{y)ds(y) 


as  m  — »  00. 

Proof:  We  introduce  regularizing  functions 
u(x)  =  j^<i>{x,y)p{y)ds{y) 

v,{x)  =  /  <^(®,y)p(y)</s{y)  +  /  <^e(x,y)p(y)d.s(y) 

Jv\B,  JrnR, 

where  (cf.  (2.2)) 

M^^y)  =  =  ^  (‘  ~2iogc  -  +  Mr) 


It  can  be  easily  verified  that  the  function 


'I»j(x,y)  =  'l',(r)  = 


4>(r)  ,  £  <  r 

<f>c(r)  ,  0  <  r  <  t 


(3.14) 


(.3.15) 


(3.16) 


(3.r 


as  a  function  of  x,y  is  a  C’  function  in 
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Assuming  that  |;^(y)l  <  M  we  have  now,  for  e  <  6, 

\v{x)  -  v,{x)\  <  ^  ^i^2log£  -  |p(y)|<^5(y) 

M  f 

<  —  /  2(logt  -  log?-)  +  -  1  My) 

4?r  JrnBt 

^  (3.18) 

2iog^  +  ^-i  v^i  +  {ny)fdu 

21og-  +  ^-l  jl  +  r“|rf^i  <  Celoge -+ 0 

r  e 

as  £  — ?  0.  Then  v  is  the  limit  of  the  uniformly  convergent  sequence  of  continuous  functions 

Vc  and  therefore  is  continuous.  B 

Lemma  3.  Let  us  suppose  all  assumptions  of  Lemma  2  hold  and,  moreover,  the  Holder 
condition  on  p  over  F  is  fulfilled.  We  have  then 

J^-^iXm,y)piy)ds{y)  — ►  ^^-(a;,y)p(y)(/.s(y)  +  ^p(x)  (3.19) 

1/  1/ 


Proof: 


_^^^(®m,y)p(y)<^s(y)  =  J^-^{xm^y)\p{y)  -  p{x)]ds{y)  +  p{x)  J^-^{xm,y)ds{y) 

=  “  p{^)]ds{y)  +  p{x)  Ay(f>{Xm,y)ds{y)  {[ 

=  /  4^{^m,y)\p{y)  -  p{x)]ds{y)  -  p{x)k^  f  <^(x,„,y)ds(y) 

Jr  any  Ju 


Next  we  take  functions  similar  to  those  in  Lemma  2; 


^(®)  =  I  ^{x^y)[p(y)  -  p{x)]ds{y 

Jr  any 

4(®)  =  f  T. —  (a:,y)(p(y)  -  p(a:)]d^i(y) 

Jr\B.  any 

c\  / 

+  /  Tj— ^(a:.y)[/?(y)  -  p(x)]t/s(y) 

JrnB,  any 


(3.20) 


(3.21) 
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where  (^e(x,y)  is  defined  by  (3.16). 


m®)  =  [ 

JrnB, 


=  Lb.  2^r  +  *»('■)  + 


27re2 


ln(y) -rl 


|p(y)  -p(®)|ds(y) 


|p(y)  -p(®)|ds(y) 


(3.22) 


<-r\--- 

2%  j-e  I  r 


|cos^(y)|r"<i^i  <  Ce”" 


where  i9(y)  is  the  angle  between  n(y)  and  r  (see  Fig.  4)  and  C  denotes  different  constants. 
Then  u  is  the  uniform  limit  of  and  is  continuous.  We  have  obtained  then 

y)[piy)  -  p{x)]ds{y)  — ^  y)[p(2/)  -  p(®)]‘^-s(y)  (3.23) 

2/  y 

In  turn,  by  Lemma  1, 


mi  y)dy  — »  k‘^p{x)  J^(i>{x,y)dy 
=  -P(®)  /  y)ds(y)  -  ^p(x) 


(3.24) 


Finally,  (3.20),  (3.23),  and  (3.24)  imply  (3.19).  I 

Remark.  The  proof  of  inequality  (3.22)  does  not  depend  upon  the  assumption  that 
®  €  r.  It  does  not  use  even  the  fact  that  |cos/3(y)|  <  Cr“  for  x  E  F  (Proposition  3).  For 
this  reason  estimation  (3.22)  does  not  depend  upon  x  and  convergence  of  Wj  is  uniform. 

Proposition  3.  If  the  first  derivative  of  the  function  /  which  describes  the  boundary  F 
(3.13)  fulfills  in  a  given  interval  (^1,^2)  the  Holder  condition  with  exponent  a,  then  there  is 
a  constant  C  >  0  such  that 


I  cos^(®)|  <  C7r",  I  cosy0(y)j  <  Cr“,  |nj.  -  n^l  <  Cr°  (3.25) 


in  this  interval. 

Proof:  We  establish  the  tangential-normal  coordinate  system  as  in  Lemma  2  (Fig.  5). 
We  have  the  following  relations 


I  cos  /3(®)| 


h  \m{y))\ 

r  r 


s‘(y)/'({)| 

r 


< 


(3.26) 


The  proof  for  cos/9(y)  is  analogous. 


of  ii{x)  upon  r. 


|3(y) 


The  vector  d  =  is  presented  in  Fig.  5. 


|dl  =  2 


.  7 
sm- 


=  2  sin  [^(x)  +  I3{y)  -  ;r]^  j  =  2  |cos  [0{x)  +  /3(y)]) 


=  \/2^cos  (/3(x)  +  /3(j/))  +  1  =  \/2\J cos ^{x) cos ^{y)  —  sin/3(x)  sin  p{y)  +  1  <  Cr° 


because  of  (3.26). 

Lemma  4.  By  the  assumptions  of  Lemma  3 


B 


^  (®m,  y)  p{y)My)  — »  ^  y)piy)My)  -  (3-27) 


t®  unx 

where  is  the  outward  normal  vector  at  x  e  F. 

Proof; 

=  /r 

d(t> 


d<j>  .  .  d<j>  .  . 

"a  i^miy)  d"  (®fn»y) 

071®  071  y 


p(y)<^5(y) 


(3.28) 


/any**”’  y)p(y)ds(y) 


The  definition  of  <f>  implies 

d4> 


dv  dv  dd> 


from  which  follows 


1  =  1 


=  -”.(y)]  (3.29) 

i=l 


Quite  analogously  to  the  proof  of  (3.23),  with  the  use  of  (3.29)  and  Proposition  3,  the 
convergence 


I 


d<f> 

dn 


X 


r  X  94>  , 

(®miy)  "h  (®miy) 


y 


p(y)<^5(y)  -»  ^ 


d(f> 


d4> 


an. 


p(y)‘^^(y) 


ma>  be  established.  This  with  Lemma  3  gives  (3.27). 
Lemma  5.  By  the  assumptions  of  Lemma  3 


(3.30) 

B 


-f  -^{Xm,y)p{y)ds{y)  — ^  -/  ■^{x,y)p{y)ds(y)  (3.31) 

Jr  OTx  Jr  OTx 

y)p{y)ds(y)  =  ^  y)piy)d4y)  (3-32) 


and 


d 

dTx 
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Proof:  For  the  beginning  let  us  consider  the  following  Cauchy  Principal  Value  integral: 
£  +  ^£  I  (losi)  dy  (3.33) 


(3.35) 


(3.36) 


(3.37) 


£  £  ('“« (;)) = ii-”o  £  ('°« (;))  ''»+££  ('°*  ?) 

(3.34) 

=  !L"3  [-  7  -  £  7  =  !i3  [/  ‘  -“T  +  i“  t]  =  “ 

(3.33)  and  (3.34)  imply 

The  same  result  can  be  obtained  when  we  replace  <!>  by  (cf.  (3.17)): 

We  define  once  more  the  “regularizing”  functions 

*  (3.37) 

^  o  lp(®.y)p(yM5(y) 

yr\fl,  ar®  JrnB,  otx 

where  rx  is  a  tangent  vector  at  x  (fixed),  and  assume  that  boundary  F  is  defined  by  function 
/  (cf.  (3.13)).  Then 

\t(x)-t^{x)\=  I  \^(x,y)-^{x,y)  p{y)ds{y) 

JVnB,  \(jTx  OTx 

-  !  b(y)-p(®)]<^^(y) 

yrnB,  {OTx  OTx  J 

+  /  p(®)<^^(y) 

JVnBt  I  OTa;  OTx 

(3.: 

~  Is  v^i  +(/'(y))^b(j/)-p(®)]‘^y 

+ £*  { [1^  •*’"> "  ®>]  -  [1^  <*’»*  ■ 

C  f  _  1 


(3.38) 
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The  last  integral  is  equal  to  zero,  because  of  (3.35)  and  (3.36),  and  the  remaining  terms  are 


<  2  r  1;^  -  ^  (1  +  Cr“)  Cr“dr  +  2M  f  1;^  -  — 
Jo  \2'Ke^  2irr  '  Jo  \2ire^  2;rr 


(I  +  rr"  -  l)dr  <  C£“ 

(3.39) 

This  inequality  proves  the  uniform  convergence  of  to  <,  then  t  is  a  continuous  function  and 
(3.31)  holds. 

In  the  next  step  we  shall  prove  the  existence  of  the  Cauchy  Principal  Value  Integral  on 
the  right  side  of  (3.32).  Let  us  denote  by  5  a  tangent  line  at  ®(^i  —  axis)  and  by 
Bgy — two  balls  with  radii  £\,£2i  centered  at  x. 

According  to  (3.35), 

Jsn{B.^\B.j)  dra 


{x,y)pix)dy  =  0 


By  definition  of  the  Cauchy  Principal  Value  Integral, 

-f-^i^,y)p{y)d4y)  =  I  ^ix.y)piy)My) 

Jr  OTx  Jr\B^^  OTx 

f 

+  .  .^{^,y)piy)My) 


(3.40) 


(3.41) 


\lr  (  X  =  1/ /  ,  .^{x,y)\piy)  - pi^)]My) 

\Jrn(B..\B,^)  OTx  I  \J^^\(B.^\B.J)  OTx 


/rn(B,,\B«j) 


+ 


d<f> 

)  dTx 

d(t> 
)  Otx 


I  .  .^{^^y)p{^)ds{y)  -  / 

•/rn(B,j\B*2  J  UTx  Js 

lsn[B,^\Bt 


d(i> 

hx 

d(f> 

5n(B,,\B,j)  Otx 


(x,y)pix)dy 


{x,y)p{x)dy 


<  I 

Jsn{B.^ 


1 


\B.j)  27rr 


dr 


dTj 


Cr^y/^+{f'(y)fdy  (3-12) 


<C(£?-£f) 

(3.42)  proves  the  existence  of  the  integral  (3.41)  in  the  Cauchy  Principal  Value  sense. 
Functions  and  (cf.  (3.37)  and  (3.15))  are  bounded  and  continuous,  moreover 

t,(x)  =  ^rv(x)  (3.43) 

()Tx 


If) 


We  have  proven  in  (3.39)  that  they  are  uniformly  convergent  to  t(x)  and  gf-v(x)  V  x  €M 


Then 


and  (3.32)  holds. 


i{x)  =  - — u(x)  VxGJE^ 

C/Tx 


4  Hypersingular  Formulation 

Let  us  begin  with  (2.12).  If  we  apply  to  it  Proposition  2,  we  obtain 

f  d<t>  ,  \  /  \j  f  \  \ 


(4.1) 


Proposition  1  implies 

dp  ^  f  d<f> 


^■(aJm)  =  J^-^iXm,y)-^iy)ds{y)  +  k^  J^TxTy<i>{x„,y)p{y)dsiy) 


dp 


dp" 


(4.2) 


Now  we  may  pass  to  the  limit  with  x..,  -♦  a:  6  F,  deriving  the  advantage  of  Lemmas  2-5: 
9P  f  y_  dp 


drii 


(«)  =  J^-^{x,y)-^{y)ds{y)  +  P  J^TxTy<i>{x,y)p{y)ds{y) 


ry 

dp 


1  dp 


dp" 


(4.3) 


The  product  differentiation  rule  implies 


_d_ 

drx 


<i>{x,y)^{y) 


dp 

^y 


^(x,s)^(»)+4.(x,»)g^(p)  =  ^(>=.»)^(»)  (4.4) 


and  this  allows  (4.3)  to  convert  into 
\  dp  .  .  f  d 


2  dn 


x'**  =  JrOTi 


<X(®,y)  J^(y) 

OTy 

dp 


ds{y)  +  k-^  j^rxry<l>(x,y)p{y)ds{y) 


fdd  \  dp  ,  ,  . 

-  r: —  ®sy)Tj-^(y)d25(y)  +  i 

Jr  Orix  ony  Orix 
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and  making  use  of  Lemma  5 


1  dp 

2d^ 


(®) 


d  t  d 

jr  I  +  k'^  I  rxTy<j>{x,  y)p{y)ds{y) 

Sf 

/r  ^  e  '' 


(4.5) 


This  is  the  final  form  of  the  hypersingular  formulation  of  (1.1)-(1.3). 

To  obtain  the  variational  formulation  of  (4.5)  let  us  multiply  both  sides  by  a  test  function 
q  and  integrate  over  F,  using  simultaneously  (3.8): 


1 

2 


dp  ,  .  dq 


Ir  Jr  y)piy)Q{x)ds(y)M^) 

-  Ur 

r  dp'^'^ 

+  Jr 

for  all  admissible  test  functions  q. 

If  we  take  a  linear  combination  of  (2.5)  and  (4.5)  with  real  and  imaginary  coefficients,  we 
obtain  the  Burton-Miller  formulation  of  the  Helmholtz  exterior  boundary-value  problem: 


2^^*^  +  -IF-  ^ 


(X)  = 


(x,y)p(y)  -  <P{x,y)-^{y) 


dn 


y 


dn 


y 


ds{y) 


+ 


(1 


k 


—  I  ^  ^  y)-^{y)My)  +  rxTy4>{x,y)p{y)ds{y)  (4.7) 

\  y 


for  Q  G  [0, 1]. 

Multiplying  both  sides  of  Eq.  (4.7)  by  a  sufficiently  smooth  test  function  q{x)  defined 
on  the  boundary  F  and  integrating  it  over  F  will  result  in  a  variational  formulation: 


IS 


Find  p  and  €  ^(r)  such  that 


®) 


(1- 


\  j^p{x)q{x)dx  +  |^^(a:)9(x)dx 

=  “/r/r 

~  {"  /r  /r 

- /r /r 


ds(y)ds(x) 


and 


dp 


(4.8) 


for  all  q  G  //(F)  where  //(F)  is  a  Hilbert  space  on  F. 

If  the  boundary  F  is  smooth  then  it  is  also  known  from  the  weak  form  of  the  boundary 
value  problem  that  the  appropriate  space  for  Helmholtz  differential  equation  is  that  p  € 
//2(F)  and  1^  G  //"2(r),  but  for  this  special  problem,  we  require  that  p  and  are  in  the 
same  space.  Thus,  we  here  assume  that  H{V)  —  //2(F)  in  the  variational  formulation  (4.8). 


5  Numerical  Implementation 

In  what  follows,  only  the  Burton-Miller  approach  is  considered.  The  claissical  integral  formu¬ 
lation  and  the  hypersingular  integral  formulation  can  be  obtained  by  simply  selecting  q  =  1 
and  0  =  0. 


Galerkin  Approximation 

Approximating  p,  and  q  in  (4.9)  with  linear  combinations  of  the  same  basis  functions 

Pi^)  =  H  P^Ck(x),  ^(x)  =  efc(x),  q{x)  =  ^  f/efc(x),  (5.1) 
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where  Ck  =  €*(*)  are  real-valued  basis  functions  and  p*',  and  q'‘  are  complex  coeffi¬ 
cients,  leads  to  a  system  of  (complex)  linear  equations  for  p*  and  .  e  =  constant  allows 

us  to  eliminate  by  setting  =  ep*  from  the  calculations  and  leads  to  a  final  system 

of  equations  of  the  form 

N 

^afc,p*=  =  6i,  (5.2) 

fc=i 

where  the  load  vector  is  given  by 

6,  =  I  (5-3) 


and  the  stiffness  matrix  is  defined  by 


{1  /■  f  f 

^  J^ek{y)ei{x)ds(x)  -  J^ekiy)-^^^{x,y)ds(y)ei{x)ds{x) 

+£  /  f  ek{y)^{x,y)ds{y)ei{x)ds{x)\ 

\  ^  (5.4) 

x{x,y)ds{y)ei{x)dsix)  +  €  ^^efc(y)^|^(x,y)ds(y)e,(x)ds(*)}  . 
h~p  Approximation 

A  one-dimensional  version  of  the  h-p  approximation  discussed  in  [4]  has  been  implemented. 
The  principal  features  of  the  method  are  as  follows: 

•  each  element  has  at  most  three  nodes,  two  endpoints  starting  with  linear  approximation 
and  a  middle  node  for  orders  of  approximation  p  >  2  (this  spectral  order  p  should  not 
be  confused  with  the  pressure), 

•  hierarchical  shape  functions  discussed  in  [7]  are  implemented, 

•  the  approximation  is  continuous.  The  linear  degrees  of  freedom  are  common  for  neigh¬ 
boring  elements,  and 

•  instead  of  quadratic  elements  used  in  the  previous  research  (see  [3]),  cubic  spline  el¬ 
ements  are  used  for  the  approximation  of  geometry.  As  shown  in  Fig.  6.  the  cubic 
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Figure  6;  Cubic  spline  interpolation. 


spline  interpolation  is  defined  by  the  map: 

+  X2<f>2iO  +  +  <*2^2(0 

3/(0  =  3/i<^i(0  +  y2<f>2U)  +  i'i^i(0  +  ^2 '^2(0 


(5.5) 


where  x(^)  and  </(  j  are  curvilinear  coordinates  with  the  master  coordinate  ^  €  [—1,1], 
and  ^Jr  jy2)  are  the  coordinates  of  the  two  end  points,  and  'I'r,  and  'I'j 


ibic  shape  functions  of  cubic 

spline,  i.e.. 

MO  = 

K  +  l)'(2-{) 
4 

MO  = 

«■(«)  = 

4 

«,({)  = 

(e-l)'(2  +  0 


4 

\2  , 


(e  +  i)0e-i) 


(5.6) 


The  map  is  C'-continuous  whenever  ai,<22,^>i  and  62  in  (5.5)  satisfy 


/  j 

yi  =  T-  and 

0l 


,  «2 


(5.7) 


where  y\  and  j/j  represent  tangential  derivative  ^  at  end  points  1  and  2,  respectively.  The 
curve  defined  in  (5.5)  will  pass  through  the  center  point  (13,3/3)  at  ^  =  0,  if  01,02,61  and  62 
in  (5.6)  further  satisfy 

1111 

(5.8) 


1111 

-Xi  +  -X2  +  -Oi  —  -02  —  X3 

1  1  I,  1, 

TiUl  +  0^2  +  761  —  762  —  7/3 

2  2  4  4 


(5.9) 
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(5.10) 


Now,  ai,a2,6i  and  62  can  be  determined  by  solving  (5.7)-(5.9)  and  we  have 

ai  =  - —  [(4^3  -  2yi  -  2^/2)  -  ^2  (4^3  -  2i,  -  2x2)) 

i/l  ~  i/2 

02  =  — -7  [(4?/3  -  2i/i  -  2j/2)  -  y'l  (4x3  -  2xi  -  2x2)]  (5-11) 

Vi  ~  i/2 

bi  =  77^  1(42/3  -  2j/,  -  2y2)  -  (4x3  -  2x,  -  2x2)]  (5.12) 

Vi  ~  Vi 

b2  =  -T^  [(4y3  -  2yi  -  2^/2)  -  y\  (4x3  -  2x,  -  2x2)]  (5.13) 

2/1  ~  2/2 

The  cubic  spline  defined  in  (5.5)  with  ai,a2,bi  and  62  determined  by  (5.10)-(5.13)  can 
fit  three  points  at  each  element  exactly  and  provide  C*-continuity  over  the  whole  boundary. 

Singular  Integration 

The  only  singular  integral  involved  in  (4.9)  is  a  weakly  singular  (logarithmic)  integral.  Let  Ki 
and  A'2  be  two  elements.  A  typical  contribution  to  the  element  stiffness  matrix  corresponding 
to  A'l  and  A'2  is  of  the  form 

Xii^Wx)  (5.14) 

where  xi  shape  functions  for  element  Kj,j  =  1,2.  Only  when  A'l  =  A'2,  the  kernel 
function  <f){x,y)  is  singular.  Thus,  in  case  of  A'l  ^  A'2,  both  outer  and  inner  integrals  are 
evaluated  using  Gaussian  integration  with  the  number  of  integration  points  iVj  =  p  +  1.  In 
the  case  A'j  =  A'2,  the  outer  integral  is  again  evaluated  using  Gaussian  quadrature  with 
N\  =  p  +  1  integration  points.  Integration  of  the  inner  integral  with  logatrithmic  singular 
type  of  (t>{x,y)  in  (5.14)  is  illustrated  in  Fig.  7.  For  every  Gaussian  integration  point  used  to 
evaluate  the  outer  integral,  the  element  is  divided  into  two  parts:  integrate  from  one  endpoint 
to  the  singular  point  and  from  the  singularity  to  the  other  endpoint  and  the  integration  is 
performed  separately  on  each  part.  Each  of  the  two  subintegrals  is  converted  into  an  integral 
from  0  to  =  1,2,  with  respect  to  the  distance  between  an  integration  point  and  the 
singular  point.  A  special  integration  rule  for  logarithmic  kernels  of  Gaussian  type  is  then 
used  with  the  number  of  integration  points  iV,  =  p  +  1.  An  additional  complication  arises 
due  to  the  lack  of  an  explicit  relationship  between  the  integration  variable  (distance  ,s)  and 
the  master  element  coordinate,  say  y.  This  solved  in  practice  by  starting  with  an  initial 
value  for  tj  corresponding  to  a  linear  element  and  performing  Newton- Raphson  iterations. 
In  practical  calculations,  the  number  of  iterations  seldom  exceeds  3. 


6  Numerical  Experiment  1 


The  classical  scattering  problem  (see  Fig.  8)  for  a  plane  wave  impinging  on  a  massless 
cylindrical  body  with  radically  elastic  response,  j.e., 

dp 

^  =  ep  (0  <  £  <  oo:  a  constant)  on  F  (6.1) 

is  used  to  test  the  method  and  code.  Using  the  polar  coordinates  (r,  0),  the  exact  solution 
to  the  test  problem  is  of  the  form 

p  =  p*""  +  p*  (6.2) 

where  the  incident  pressure 

ptnc  _  p.^^exp{ikr  cos6)  =  P,„cexp(itzi)  (6.3) 

with  a  constant  P,„ci  and  the  scttered  pressure  p*  is  evaluated  by  the  series 

P"{r,q)  =  B,nfImikr)cos{m0)  (6.4) 


with  /fm-Hankel  functions  of  the  first  kind  of  order  m,  and  coefficients  8,^  are  given  by 

dJ^{ka]  €  ,  , 

dlkr) 

5m  =  . : -  (6.5) 

dHm{ka)  £ 

-KhT  ~  it  ’ 

where  Jm  is  Bessel  function  of  the  first  kind  of  order  m,  a  is  the  radius  of  the  cylinder,  and 


6o  =  1  and  =  2  (m  =  1,2,3, . . .)  (6.6) 

k  is  the  wave  number.  In  the  case  £  =  0  the  problem  reduces  to  the  classical  rigid  scattering 
problem. 

In  practice,  we  can  only  calculate  a  finite  number  of  terms  of  the  form  (6.5)  in  (6.4)  to 
obtain  a  solution  for  the  test  problem.  The  calculation  of  the  Hankel  and  Bessel  function 
also  introduces  some  errors.  Based  on  numerical  experience,  a  reasonable  estimation  of  the 
L^-error  of  the  solution  calculated  by  using  (6.1)-(6.5)  is  between  10“^  to  10“®. 


Forbidden  Frequencies  of  the  Test  Problem 

The  solutions  of  the  Helmholtz  formulation  and  the  hypersingular  formulation  for  exterior 
acoustic  wave  problems  is  nonunique  at  certain  forbidden  frequencies  or  wavenumbers.  The 


forbidden  wavenumbers  are  found  to  coincide  with  the  “resonant”  wavenumbers  for  a  related 
interior  problem.  In  the  two-dimensional  case  of  a  circle  of  radius  a  (the  test  problem),  for 
example,  the  forbidden  wavenumbers  of  the  Helmholtz  formulation  consist  of  values  of  k  such 
that  Jn{ka)  =  0  and  the  forbidden  wavenumber^s  of  the  hypersingular  formulation  consist 
of  values  of  k  such  that  J'^{ka)  =  0,  where  is  a  Bessel  function  of  the  first  kind  and 

order  n  and  n  is  integer. 

Example  1.  Comparison  of  solutions  for  two  kinds  of  hypersingular  formulations. 

In  this  example,  we  compare  results  solved  by  the  Helmholtz  formulation  (2.5),  the 
current  hypersingular  formulation  (4.6),  and  the  classical  hypersingular  formulation  (2.7). 
The  test  problem  is  solved  using  uniform  meshes  of  8,  16,  32,  and  64  elements  with  order  of 
approximation  p  =  1,2, ...  ,6.  The  following  data  are  chosen  in  the  test  problem: 

£  =  0  (rigid  scattering) 

P.nc  =  1 
a  =  1 
^:  =  4 

The  comparisons  of  the  T^-errors  of  the  test  problem  solved  by  three  different  formula¬ 
tions  for  p  =  2  and  p  =  3  are  presented  in  Figs.  9  and  10.  The  results  for  other  orders  of 
approximation  are  similar  and  are  not  illustrated.  Figures  9  and  10  indicate  that  the  hy¬ 
persingular  formulation  (4.6)  can  give  a  much  better  solution  than  the  classical  formulation 
(2.7).  In  Figs.  9  and  10,  the  solution  of  the  Helmholtz  formulation  is  slightly  better  than 
the  solution  of  the  hypersingular  formulation  (4.6),  but  there  is  no  substantial  difference 
between  the  solutions  solved  by  two  formulations. 

It  seems  that  the  classical  hypersingular  formulation  (2.7)  results  are  not  completely  sat¬ 
isfactory  in  practical  computation.  In  most  cases,  the  accuracy  is  much  lower  than  that  of 
the  Helmholtz  formulation.  The  low  accuracy  in  the  hypersingular  formulation  deteriorates 
the  corresponding  solution  of  the  Burton-Miller  formulation.  A  major  purpose  of  this  re¬ 
search  is  to  improve  the  accuracy  of  hypersingular  formulation  solution  and  we  have  arrived 
at  this  purpose. 

Example  2.  Illustration  of  the  forbidden  frequencies  for  the  rigid  scattering  problem. 

The  purpose  of  this  example  is  to  show  the  effects  of  forbidden  frequencies  on  formula¬ 
tions.  Theoretically,  as  long  as  the  wave  number  k  does  not  coincide  exactly  with  a  forbidden 
frequency,  the  corresponding  stiffness  matrix  is  nonsingular  and  the  approximate  problem 
has  a  unique  solution.  In  practice,  for  wave  numbers  k  close  to  the  forbidden  frequenrics. 
the  condition  of  the  matrix  deteriorates  and  the  quality  of  the  solution  drastically  decreases. 
The  situation  is  illustrated  in  Figs.  11  and  12  showing  the  variation  of  the  approximate  and 
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exact  solutions  at  field  point  (—3,0)  as  t  function  of  wave  number  k.  The  following  data 
have  been  chosen: 


e  =  0  (rigid  scattering) 

Pine  =  1 
a  =  1 

Afc  =  0.001 

^min  =  1,  and  A:n,ax  =  9 

For  each  k,  the  test  problem  was  solved  by  the  Helmholtz  formulation  (shown  in  Fig.  11) 
and  the  hypersingular  formulation  (shown  in  Fig.  12)  on  the  same  uniform  mesh  of  8 
quadratic  elements.  Both  formulations  lose  stabilities  at  certain  frequencies.  Figure  13  shows 
calculations  of  the  Burton-Miller  approach  (based  on  the  current  hypersingular  formulation 
(4.6)).  The  formulation  delivers  consistently  stable  results  for  the  whole  range  of  wave 
numbers  k. 

The  solution  produced  by  the  current  hypersingular  formulation  (4.6)  proves  very  stable. 
With  Ak  =  0.001,  the  numerical  solution  in  Fig.  11  only  loses  stability  at  3  forbidden 
frequencies.  With  A^;  =  0.002,  kj^^  =  1,  and  =  8,  the  same  experiment  in  [3]  showed 
that  the  numerical  solution  based  on  (2.7)  lost  stability  at  9  forbidden  frequencies  (see  Fig. 
14).  This  simply  implies  that  the  condition  of  the  stiffness  matrix  formed  by  the  current 
hypersingular  formulation  is  better  than  the  condition  of  the  matrix  formed  by  the  cleissical 
formulation  when  near  a  forbidden  frequency. 

7  Convergent  Rate  and  A  Posteriori  Error  Estimation 


Convergence  of  the  Galerkin  method  for  boundary  integral  equations  has  been  carefully 
analyzed  in  a  series  of  works  by  Wendland  [10]  and  p-  and  /ip- versions  of  boundary  elements 
are  discussed  by  Rank  [7]  and  Stephan  and  Suri  [8,9]  and  a  brief  discussion  was  given  in 
[6].  Here,  we  follow  the  same  approach  in  [3]  together  with  a  posteriori  estimates  and  use 
Z^-norm  of  error  to  measure  convergent  rates  of  the  formulations,  (2.5),  (4.6),  and  (4.7). 

The  error  estimates  used  in  this  research  are  based  on  the  L^-norm  residual.  For  any 
approximation  solution  p/,,  we  define 


r(x)  =  ph{x)  -  2^ 


p(y) 


d<t>{x,y)  dpkiy) 


dn^ 


dn^ 


<t>{x,y) 


,ls  -  2p'^^{x) 


(7.1) 


as  the  residual  for  the  Helmholtz  formulation  and  then  the  L^-norm  of  residual  is  calculated 
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by 


=  J^(r{s))'^ds  =  J2j^  (r(5)) 


ds 


N  Ge 

EE' 

e=l 1=1 


K) 


(7.2) 


where  N  is  the  number  of  elements,  Ge  is  the  aumber  of  Gauss  points,  wl  represents  the 
weight  of  the  Gauss  integral,  and  s®  represent  the  coordinate  of  the  Gauss  points.  Wo  use 
Ge  =  p  +  1  in  (7.2)  for  each  element,  p  being  the  order  of  approximation  in  each  element. 


Based  on  the  analysis  in  [3],  we  have  the  following  theorem  for  the  residual  defined  in 
(5.1). 

Theorem: 


Let  e  =  p  —  ph  denote  the  approximation  error  inherent  in  p)^.  Then,  for  the  rigid 
scattering  problem,  if  k  is  not  a  forbidden  frequency  of  (2.5),  there  exist  positive  constants 
C\  and  Ci  such  that  for  the  residual  r  defined  in  (7.1)  and  error  e  there  holds 


Ci||rli<lle||<C2llrl| 


(7.3) 


where  ||  •  ||  represents  the  L^-norm  on  F  and  p  is  the  exact  solution.  In  (7.3),  e  can  be  the 
error  from  the  solution  of  the  Helmholtz  formulation,  the  hypersingular  formulation,  or  the 
Burton-Miller  formulation.  Furthermore,  if  e  is  the  error  from  the  solution  of  the  Helmholtz 
formulation,  then 


lim 

h—Op—oo 


where  h  is  the  maximum  size  of  the  mesh  and  p  is  the  order  of  approximation. 


(7.4) 


In  case  t  is  a  forbidden  frequency,  there  still  exists  a  positive 


constant  C  such  that 


IlHell  <  Cllrll  (7.5) 

where  H  is  the  L^-orthogonal  projection  of  p  onto  the  orthogonal  complement  of  N{A),  A  is 
defined  as  the  integral  operator  A 

Ap  =  ^p{x)  - 

and  N{A)  is  the  null  space  of  A.  3 

Compared  with  other  forms  of  error  estimates  proposed  in  the  boundary  element  method 
(see,  e.g.  [6,7]),  the  L^-norm  of  the  residual  derived  here  is  simple  to  calculate  and  the 
computed  results  agree  with  the  theory  quite  well  in  [3]. 


8  Numerical  Experiment  2 

Example  3.  Rates  of  convergence  for  the  Helmholtz  formulation. 
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Selecting 


e  =  0  (rigid  scattering) 

P  =1 
a  =  1 
it=  15 

The  test  problem  was  solved  using  uniform  meshes  of  8,  16,  32,  and  64  elements  with  order 
of  approximation  p  =  1,2, ...,6.  The  corresponding  rates  of  convergence  for  /i-refinements 
at  fc  =  15  are  summarized  in  Figs.  15  and  16  and  for  /^refinement  in  Fig.  17.  The  rates 
of  convergence  for  the  A-method  are  very  close  to  the  theoretical  rate  p  +  1,  increasing  in 
particular  with  the  order  of  approximation  p.  The  rates  of  convergence  for  the  p-extensions 
seem  to  deliver  the  superlinear  (exponential)  rate  of  convergence. 

Example  4.  Rates  of  convergence  for  the  hypersingular  formulation. 

With  the  same  data  and  the  same  h-p  meshes  as  in  the  first  example,  the  problem  was 
solved  once  again  using  the  hypersingular  formulation  given  in  (4.6).  The  results,  summa¬ 
rized  in  Figs.  18,  19,  and  20,  arc  slightly  less  satisfactory  than  those  in  the  first  example. 
The  observed  rates  of  convergence  for  the  /i-method  are  almost  not  changed  for  p  >  3.  The 
p-method  shows  no  signs  of  exponential  convergence.  The  reason  for  these  suboptimal  rates 
of  convergence  is  still  not  clear  at  this  writing,  supposedly  because  of  the  domination  of  the 
error  due  to  the  geometry  approximation.  The  results  produced  by  the  hypersingular  formu¬ 
lation  are  more  sensitive  to  the  geometry  approximation  than  those  obtained  by  using  the 
Helmholtz  formulation;  thus,  the  cubic  spline  interpolation  in  geometry  may  not  be  accurate 
enough  for  p  >  3. 

Example  5.  Global  effectivity  of  the  L^-residual  estimate. 

We  begin  with  a  simple  comparison  of  the  L*-error  and  T^-residual  for  the  test  problem 
with  data  and  meshes  discussed  in  Examples  1  and  2.  The  effectivity  of  the  L^-residual 
estimate  is  measured  by  an  effectivity  index  which  is  defined  by 


effectivity  index  = 


M 

Hell 


(8.1) 


Here,  ||  •  ||  is  the  norm  of  7^^(r),  r  is  the  residual  defined  in  (7.1),  and  e  =  p^  —  p/,  where 
Ph  is  the  numerical  solution  of  the  test  problem  and  pg  is  the  solution  solved  by  (6.1)-(6.6). 
Figure  21  summarizes  results  for  the  Helmholtz  formulation.  Fig.  22  for  the  hypersingular 
formulation,  and  Fig.  23  for  the  Burton-Miller  formulation.  The  following  conclusions  can 
be  drawn. 


•  For  the  Helmholtz  formulation  the  L^-rcsidual  seems  to  give  an  asymptotically  exact 
error  estimate  for  /(-refinements,  as  predicted  by  the  Theorem  in  (7.4). 
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•  For  large  p  and  small  size  elements,  the  Helmholtz  formulation  effectivity  index  de¬ 
creases,  supposedly  because  of  the  domination  of  the  error  due  to  the  geometry  ap¬ 
proximation  or  due  to  the  error  in  calculating  the  exact  solution  by  (6.1)-(6.6). 

•  For  the  hypersingular  formulation,  the  global  effectivity  index  stays  bounded  (between 
1  and  2  for  the  test  problem)  but,  as  expected,  does  not  converge  to  1. 

•  The  effectivity  index  of  the  Burton-Miller  formulation  also  stays  bounded  and  is  very 
close  to  1. 

Example  6.  Global  effectivity  of  the  L^-residual  estimate  for  a  forbidden  frequency. 

Figure  24  presents  the  results  of  the  error  estimation  for  the  test  problem  at  a  forbid¬ 
den  frequency  of  Helmholtz  formulation  {k  =  5.1356)  obtained  by  using  the  solution  of 
the  Burton-Miller  formulation.  The  displayed  values  of  the  global  effectivity  indices  stay 
consistently  in  the  same  range  as  in  Fig.  23  (close  to  1),  indicating  that  the  T^-residual 
corresponding  to  the  Helmholtz  formulation  can  be  used  for  the  entire  range  of  frequencies, 
provided  the  iV(/4)-component  of  the  error  is  eliminated  by  an  appropriate  formulation  (in 
this  case,  the  Burton-Miller  formulation). 

For  large  p  and  small  size  elements,  the  decrease  of  the  effectivity  index  is  observed  in  Fig. 
24.  This  is  believed  to  be  the  domination  of  the  error  due  to  the  geometry  approximation 
or  due  to  the  error  in  calculating  the  exact  solution  by  (6.1)-(6.6). 

Example  7.  Local  effectivity  of  the  T^-residual  estimate. 

The  efficiency  of  the  adaptive  method  is  directly  determined  by  the  local  property  of 
the  error  estimate.  With  the  same  data  as  in  Examples  5  and  6,  the  local  L^-error  and 
Z/^-residual  is  also  compared.  The  local  effectivity  index  in  each  element  is  very  close  to  the 
global  effectivity  index.  As  an  example.  Table  1  presents  a  comparison  of  local  Z-^-error  and 
T^-residual  is  also  compared.  The  local  effectivity  index  in  each  element  is  very  close  to  the 
global  effectivity  index.  As  an  example.  Table  1  presents  a  comparison  of  local  Z^-error  and 
Z^-residual  for  a  mesh  of  16  quadratic  elements  and  Burton-Miller  formulation  at  ^  =  15 
and  k  =  5.1356  (a  forbidden  frequency  for  Helmholtz  formulation). 
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Figure  7;  Evaluation  of  singular  integrals. 


Figure  8:  The  test  problem. 
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Figure  9:  Example  1.  Comparison  of  /-^-errors  of  solutions  among  three  formulations  for 


Figure  10:  Example  1.  Comparison  of  L^-errors  of  solutions  among  three  formulations  for 
p  =  3. 
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Figure  13:  Example  2.  Solution  at  a  field  point  as  a  function  of  wave  number  for  Bur¬ 
ton-Miller  formulation. 
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Log  N  (N:  Degree  of  Freedom) 


Figure  15:  Example  3.  Experimental  rates  of  convergence  for  uniform  /i-refincment  (p 
2.  and  3). 
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Log  N  (N:  Degree  of  Freedom) 


Figure  16:  Example  3.  Experimental  rates  of  convergence  for  uniform  /j-refmement  (p  =  4, 
5,  and  6). 
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Log  N  (N:  Degree  of  Freedom) 


Figure  17:  Example  3.  Experimental  rates  of  convergence  for  uniform  p-refinement. 


Log  N  (Degree  of  Freedom) 


Figure  18:  Example  4.  Experimental  rates  of  convergence  for  uniform  /i- refinement  (p  =  1 
2,  and  3). 
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Log  N  (N:  Degree  of  Freedom) 


Figure  19:  Example  4.  Experimental  rates  of  convergence  for  uniform  /t-refinement  (p 
5,  and  6). 
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Log  N  (N:  Degree  of  Freedom) 


Figure  20:  Example  4.  Experimental  rates  of  convergence  for  uniform  p-retinement 


Figure  21:  Example  5.  Comparison  of  global  £, ^-errors  with  L^-residual  for  the  Helmholtz 
formulation  (A:  =  15). 


N;  Degree  of  Freedom 

Figure  22:  Example  5.  Coujparison  of  global  Z/^-errors  with  L^-residuai  for  the  hypersingular 
formulation  {k  =  15). 


•42 


0  100  200  300 

N:  Degree  of  Freedom 

Figure  23:  Example  5.  Comparison  of  global  L^-errors  with  //^-residual  for  the  Burton-Miller 
formulation. 


Figure  24:  Example  5.  Comparison  of  global  L^-errors  with  Z,^-residual  for  the  Burton-Miller 
formulation  (k  =  15). 
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Table  1:  Example  6.  Comparison  of  local  Z-^-errors  and  Z^-residuals  for  a  mesh  of  32 
quadratic  elements  and  Burton-Miller  Formulation. 


Elem.  k  =  5.1356  k  =  15 
No.  Local  Effectivity  Index 


1 

0.9928 

2 

0.9991 

3 

0.9976 

1.0001 

4 

1.0017 

0.9959 

5 

0.9974 

1.0009 

6 

1.0004 

0.9970 

7 

0.9999 

0.9990 

8 

0.9968 

1.0014 

9 

1.0035 

0.9946 

10 

1.0104 

0.9959 

11 

0.9996 

1.0004 

12 

0.9905 

0.9699 

13 

1.0118 

0.9818 

14 

0.9820 

1.0067 

15 

0.9890 

0.9213 

16 

1.0138 

1.0016 

Global  Effectivity  Index 
0.9992  0.9963 
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Appendix  B 


Adaptive  Finite  Element  Methods 
For  Hyperbolic  Systems  With  Application  to 

Transient  Acoustics 


Abstract 

The  solution  of  the  hyperbolic  systems  of  equations  governing  transient  acoustics 
by  adaptive  finite  elements  and  various  finite  difference  time  discretization  schemes  is 
addressed.  Emphasis  is  placed  on  the  use  of  a  class  of  implicit  Runge-Kutta  methods 
for  temporal  approximations. 


1  Introduction 

This  paper  describes  general  high-order  adaptive  finite  element  methods  for  solving  wave 
propagation  problems  characterized  by  hyperbolic  systems  of  equations,  with  particular  em¬ 
phasis  on  problems  of  radiation  and  scattering  in  transient  acoustics.  Special  features  of  the 
present  analysis  are  numerated  as  follows: 

1.  We  accept  as  a  model  of  acoustic  wave  phenomena  the  full  hyperbolic  system  of  con¬ 
servation  laws  derived  from  a  perturbation  analysis  of  the  compressible  Navier-Stokes 
equations.  This  is  contrary  to  most  approaches  for  acoustic  modeling  which  are  based 
on  assumptions  of  a  uniform  background  flow  which  renders  these  systems  reducible 
to  a  problem  in  the  frequency  domain.  The  results  for  different  frequencies  must  then 
be  stored  and  then  superimposed  to  solve  the  fully  transient  problem.  In  practice,  this 
limits  the  application  to  cases  in  which,  at  most,  a  few  hundred  harmonics  contribute 
to  the  solution. 

The  use  of  the  full  transient  formulation  also  provides  a  basis  for  considering  the  non- 
classical  equations  of  acoustics  resulting  from  the  linearization  of  the  compressible  Euler 
equations  around  an  arbitrary  solution  to  the  incompressible  Euler  equations.  The  re¬ 
sulting  system  of  equations  is  still  linear  and  hyperbolic  but  with  variable  coefficients 
depending  on  the  back^  "umd,  incompressible  flow.  The  presence  of  space  or  time 
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varying  coefficients  renders  inapplicable  most  of  the  classical  methods  (in  particular, 
those  designed  for  treating  the  reduction  to  the  wave  equation).  Also,  by  generalizing 
the  characterization  of  the  flow  that  governed  by  the  Navier-Stokes  equations,  viscous 
effects  can  be  also  included  in  this  as  well.  Once  the  transient  formulation  is  accepted, 
additional  generalizations  such  as  these  can  be  easily  incorporated  into  an  existing 
code,  making  the  whole  approach  very  flexible. 

2.  We  assume  that  the  final  discrete  model  is  obtained  by  combining  a  time-discretization 
algorithm  with  an  h-p  FE  discretization  in  space.  In  this  sense  the  present  work  is  a 
continuation  of  our  study  on  h-p  FE  methods  presented  in  [2,8,11].  Only  linear  time 
discretization  schemes  are  considered,  but  both  the  method  of  lines  and  the  method  of 
discretization  in  time  (see  Section  3)  are  admitted  in  our  formulation. 

3.  The  study  is  directed  towards  designing  fully  automatic  self-adaptive  schemes  with 
error  control.  The  use  of  adaptive  methods  for  wave  propagation  phenomens  seems  to 
be  meaningful  only  for  the  class  of  problems  with  “short  signal”  propagation,  where 
an  initial  disturbance  moves  throughout  the  domain  but  stays  local,  and  therefore, 
requires  only  the  addition  of  extra  degrees  of  freedom  locally.  Several  fundamental 
questions  present  themselves  immediately: 

(a)  Space-Time  Consistency:  For  a  given  high-order  h-p  approximation,  what  par¬ 
ticular  time-discretization  scheme  should  be  selected?  What  are  the  desired  prop¬ 
erties  of  an  optimal  scheme,  especially  in  the  context  of  adaptive  methods? 

(b)  Error  Estimation:  Ideally,  an  error  estimate  should  estimate  errors  resulting  from 
discretization  both  in  time  and  space  variables.  How  can  this  be  implemented  into 
a  practical  analysis  strategy? 

(c)  Adaptive  Strategies:  What  adaptive  strategies  are  optimal?  Should  both  h  or 
p-refinements  be  used?  Is  there  a  need  (and  place)  for  h-p  strategies? 

The  study  of  hyperbolic  .systems  using  ideas  of  classical  spectral  theory,  applied  to 
acoustics,  provides  a  general  framework  for  addressing  these  questions  and  is  used  as 
a  basis  for  the  present  investigation. 

The  plan  of  the  presentation  is  as  follows.  Following  this  introduction,  we  summarize 
the  essential  mathematical  properties  of  the  model  in  Section  2.  Section  3  is  dev'oted  to 
a  discussion  of  two  main  discretization  concepts:  the  method  of  lines  and  the  method  of 
discretization  in  time.  A  study  on  implicit  Runge-Kutta  methods  is  presented  in  S<'ction  4 
and  some  preliminary  results  on  adaptivity  are  given  in  Section  o. 


2  Equations  of  Linear  Acoustics — A  Summary 

We  begin  with  a  brief  review  of  some  fundamental  mathematical  results  for  the  linear  acous¬ 
tics  equations.  For  a  detailed  analysis  of  the  subject,  we  refer  to  [6]. 

As  a  starting  point,  we  record  the  conservation  equations  of  isentropic,  compressible 
inviscid  flow  in  the  form  (see,  e.g.,  [5,  7]) 

P,t  +  {puk),k  =  0 

p  (ufc  +  uiUk,()  +  p,fc  =  0  ■  (2.1) 

p  =  Ap"’’,  A>0,  7>1,  l<k,  i<n  =  2  or  Z 

where  p  is  the  density,  u  =  (uk)  the  velocity  vector  (repeated  indices  are  summed),  p  the 
pressure  and  the  standard  notation  for  differentiation  is  used  (p,f  =  dp/dt,  p,k  =  dp/dxk, 
etc.). 

Linearizing  the  equations  around  the  equilibrium  state 

P  =  Po  —  const,  «  =  0  (2.2) 

and  introducing  the  sound  speed  co  defined  as 

4  =  ^  (Po)  (2.3) 

dp 
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we  arrive  at  the  classical  equations  of  linear  acoustics  in  the  form 

P,t  +  =  0 

Uk,t  +  ~  P,k  =  0 

Po 


(2.4) 


where  now  p  and  Uk  denote  perburbations  of  density  and  velocity  components  from  the 
equilibrium  state.  Introducing  the  perturbed  pressure  p  defined  by 


P  =  CoP 

we  reformulate  the  equations  in  terms  of  pressure  and  velocity 

p^t  4-  cl  pQUk,k  =  0 


Uk,t  +  —  P,k  =0 
Po 


(2.5) 


(2.G) 


It  is  important  to  note  that  the  entire  linearization  procedure  can  be  performed  around  an 
arbitrary  solution  of  the  incompressible  Euler  equations  resulting  in  the  generalized  equations 
of  linear  acoustics  (see  [7]). 

Finally,  introducing  the  nondimensional  variables 


Xk  tC(\  Uu  p 

Xk:=f.  Uk:=-f,p:=J^ 

h  L  Co  CqPq 


(2.7) 


where  L  is  a  unit  of  length,  we  arrive  at  the  nondimensional  version  of  the  equations  in  the 
form 


u^i  +  grad  p  =  0 
[  P  t  +  div  u  =0 
Introducing  the  group  variable  U  =  {u^,p)^  and  a  (formal)  operator  A 

0  grad  \  f  u  ^ 

div  0  )  \  P  / 

where  i  is  the  imaginary  unit,  we  rewrite  (2.8)  in  the  abstract  form, 

U,t  +  tAU  =  0 


(2.8) 


AU  -f 


(2.9) 


(2.10) 


Equation.s  (2.10)  are  to  be  solved  in  a  domain  Q  C  i?",  n  =  2, 3.  Typically,  two  particular 
cases  are  of  interest: 
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•  interior  problems  when  Q,  is  bounded 

•  exterior  problems  when  fi  is  a  complement  of  a  bounded  set 

In  both  cases  we  restrict  ourselves  to  two  kinds  of  boundary  conditions: 

Case  1.  Kinematic  boundary  conditions  (vibrating  boundary) 

Un^=u-n  =  Un  on  r„  (2.11) 

where  n  is  the  outward  normal  unit  to  the  boundary  and  u„  is  prescribed  velocity 
of  the  vibrating  boundary  (“solid  wall”  if  «„  =  0). 

Case  2.  Pressure  boundary  condition 

p  =  p  on  Pp  (2.12) 

where  p  is  a  prescribed  pressure  on  Pp — part  of  the  boundary  {dO,  =  Pu  U  Pp,  Pu  D 
Pp  =  0).  The  initial  boundary  value  problem  is  completed  by  specifying  initial 
condition  of  the  form 

u  =  Uq  and  p  =  Po  at  t  =  0  (2.13) 

or,  equivalently, 

U  =  Uo,  where  Uo=  («o>  Po)^  (2.14) 

In  the  case  of  homogeneous  boundary  conditions,  the  problem  can  be  cast  into 
the  Hilbert  space  formulation  as  follows. 

We  introduce 


•  The  Hilbert  space 


H  =  L^{n)  X  Z,2(fi)  =  I  JJ  Z.^(H)  1  X  L^{n) 


•  Operator  A  :  H  D  D{A)  — >  H 

where  the  domain  of  A.  D(A),  consists  of  all  vectors  U  =  (u^,p)^  such  that 


div  u  G  L^(Q)  ,  =  0  on  P„ 

pG//*(f2)  ,  p  =  0  on  Pp 


(2.15) 


Note  that  the  boundary  condition  on  p  is  satisfied  in  the  sense  of  the  tmce  theorem, 
whereas  the  boundary  condition  on  u  is  interpreted  in  the  sense  of  the  generalized 
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Green’s  formula  (comp.,  e.g.,  [13]).  For  that  reason,  the  pressure  boundary  condi¬ 
tion  is  classified  as  the  Dirichlet  boundary  condition  and  the  kinematic  boundary 
condition  as  the  Neumann  boundary  condition  for  operator  A.  Note  also  that 
the  Zi^-norm  of  the  solution  vector  U  can  be  interpreted  as  the  total  (mechanical) 
energy  of  the  field. 


Within  the  Hilbert  space  formalism,  the  initial  boundary  value  problem  can  be  reinter¬ 
preted  as  an  abstract  Cauchy  problem  for  operator  A- 


U,i  +  iAU  =  0  t  >  0 

u  =  Uo  t  =  0 

where  Uq  E  H. 

An  jFf- valued  function  of  time  U  =  U{t) 

[0,  oo)  B  t  U{t)  G  H 


(2.10) 


is  called  a  weak  solution  to  the  Cauchy  problem  if  U{t)  satisfies  two  conditions 


(i)  a  regularity  assumption 

Lr€C'([0,oo)  ;  H) 

(ii)  a  weak  form  of  the  equation  given  by 

dx  dt 

—  f  ■)dx  =  0 

Jq 

for  every  test  function 

«P  €  Co{R  ;  D{A))  nC\E-,  H) 


(2.17) 


(2.18) 


(2.19) 


Note  that  this  definition  admits,  in  particular,  solutions  in  the  d’.Alembert  sense. 

We  record  now  some  fundamental  results  concerning  operator  A  and  the  existence  and 
uniqueness  of  weak  solutions  U . 

•  Operator  A  is  self-adjoint  (in  the  complex  sense)  and  therefore  its  spectrum  lies  on 
the  real  line  and  consists  of  a  point  spectrum  (eigenvalues)  and  continuous  spectrum 
(generalized  eigenvalues)  only. 
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•  For  the  interior  problems  (f2  bounded),  the  spectrum  of  A  consists  only  of  eigenvalues 
symmetrically  located  on  two  parts  of  the  real  axis  and  “escaping"’  to  infinity  (.d  is 
unbounded) 


^{■^)  —  {0>  -^2,  —-^2)  •  •  •} 


0  <  Ai  <  A2  <  . . .  <  A„  oo 


(2.20) 


Except  for  the  0-eigenvalue,  all  eigenvalues  are  of  finite  multiplicity  and  the  eigenvectors 
corresponding  to  A„  and  — A„  are  complex  conjugate  to  each  other.  All  eigenspaces  are 
orthogonal  to  each  other. 


•  For  the  exterior  problems,  the  discrete  spectrum  reduces  to  the  single  eigenvalue  0 
while  the  rest  of  the  real  axis  constitutes  the  continuous  spectrum. 

•  In  both  cases,  the  eigenspace  corresponding  to  zero  eigenvalue,  i.e.,  the  null  space  of 
operator  A,  is  of  inifinite  dimension  and  it  contains  all  incompressible  velocity  fields. 
More  precisely,  it  has  the  form 

M{A)  =  6  H  :  div  m  =  o|  (2.21) 

where  p  is  an  arbitrary  constant  for  ft  bounded  and  Fp  =  0  or  p  =  0  otherwise. 

•  In  both  cases,  operator  A  admits  a  classical  spectral  decomposition.  In  the  case  of 
bounded  ft,  it  reduces  to  the  series  representation 

00 

AU=  (2.22) 

trs— 00 

where  dP,-  is  the  orthogonal  projection  on  the  eigenspace  corresponding  to  A,.  For 
exterior  problems, 

AU  =  r  \dPx{U)  (2.23) 

J —oo 

where  the  integral  is  understood  in  the  Riemann-Stieltjes  sense  (see,  e.g.,  [14])  and  P\ 
denotes  the  spectral  family  of  A. 


•  A  weak  solution  U  exists  and  it  is  unique.  Moreover,  it  is  of  the  form 

U{t)  =  r  e-'^^dPx{Uo)  (2.24) 

J  —  00 

where  the  integral  reduces  to  the  series  representation  for  the  interior  problems  and  it 
is  understood  in  the  Riemann-Stieltjes  sense  for  the  exterior  problems.  From  the  form 
of  the  solution,  it  follows  in  particular  that  the  energy  is  conserved 

l|e(0ll  =  IIUoll  V  ( >  0  (2.'->5) 
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•  If  the  initial  condition  function  Uq  satisfies  an  additional  regularity  assumption 


t/o  €  D{A)  (2.26) 

then  the  solution  U  €  C'^((0, oo),  Jf)nC((0,  c»),  £)( A).  We  say  then  that  U  is  a  strict 
solution  to  the  problem. 


3  Discretization 

Any  approximation  of  a  transient  problem  must  involve  discretization  both  in  time  and  space 
variables.  Although  a  simultaneous  discretization  in  all  variables  is  possible  (e.g.,  space-time 
finite  elements),  we  will  adopt  the  assumption  here  that  the  final  approximation  is  obtained 
by  using  finite  differences  in  time  and  finite  elements  in  space  variables. 

Two  different  approaches  are  possible.  In  the  classical  method  of  lines,  an  approximation 
in  space  variables  converts  the  original  initial-boundary  value  problem  into  a  system  of 
ordinary  differential  equations  (ODEs)  which  next  is  discretized  in  time  using  one  of  many 
time  integration  schemes  for  ODEs.  An  alternative  procedure  known  as  the  method  of 
discretization  in  time  (also  called  the  Rothe’s  method,  see  [12]),  consists  of  the  same  two 
steps  but  done  in  the  reverse  order.  By  discretizing  in  time  first,  the  initial-boundary  value 
problem  is  converted  into  a  sequence  of  boundary-value  (-like)  problems  which,  in  turn,  give 
a  basis  for  a  spatial  approximation  and,  consequently,  a  fully  discretized  scheme. 

Both  procedures  are  symbolically  depicted  in  Fig.  1.  If  A  denotes  the  original  operator, 
by  Afi  we  mean  its  discrete  counterpart  obtained  as  a  result  of  approximation  in  space 
performed  first  and  T{Ah)  will  denote  a  fully  discrete  transient  operator  resulting  from  the 
time  discretization  of  the  system  of  ODEs.  In  the  notation  we  have  restricted  ourselves  to  one 
time  step  methods  assuming  that  the  approximate  solution  at  time  level  is  obtained 
using  the  solution  from  the  previous  time  level,  only.  A  more  general  situation  is,  of 

course,  possible. 

In  Rothe’s  method  the  transient  operator  T  =  T{A)  is  defined  on  the  infinite- 
dimensional  level  and  only  an  appropriate  discretization  in  space,  performed  ne.\t,  converts  it 
into  a  fully  discrete  operator  Th.{A).  Surprisingly  enough,  even  if  exactly  the  same  approx¬ 
imation  in  time  and  space  variables  is  used  and  even  for  linear  equations  like  the  acoustics 
problem,  the  two  procedures  may  result  in  two  different  methods. 

We  |:)roceed  now  with  an  example  of  such  a  situation.  Assume  tliat  U  is  a  strict  solution 
to  the  problem.  Taking  (2.10)  and  multiplying  both  sides  by  the  (complex  conjugate'  of)  a 
test  function  V  e  D{A)  we  arrive  at  the  equivalent  [D[A)  is  dense  in  H)  weak  form  of 
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(2.10) 


~(U,V)+HAU,V)  =  0  VV€D{A)  (3.1) 

at 

Restricting  ourselves  now  to  a  finite  dimensional  subspace  X^^  of  D{A),  we  can  formulate 
the  corresponding  finite  dimensional  problem  as 

j^{U^,V^)  +  iiAU^,Vf,)  =  0  VV.gX^  (3.2) 

or  equivalently 

j^UH  +  iAkUH  =  0  (3.3) 

where  the  approximate  operator  Ah  '■  X/j  — is  defined  as 

Ah  =  PhoA\x^  (3.4) 

provided  Ph  denotes  the  orthogonal  projection  in  H  onto  Xh.  Note  that  the  approximation 
Ah  preserves  the  original  properties  of  A,  i.e.,  Ah  is  self-adjoint  on  the  Xh, 

{AhUh,Vh)  =  {Uh,AhVh)  V  Uh,Vh  e  Xh  (3.5) 


The  system  of  ODEs  (3.3)  has  yet  to  be  discretized  in  time.  This  can  be  done  for  instance 
by  using  the  following  finite  difference  formula  of  second  order 


U h{t  +  At)  —  a  Uh,tt{i  +  At) 

At^ 

=  U,(t)  +  AtUhAt)  +  (1  -  Uh,tt{t)  +  0(A«") 


(3.6) 


Making  use  of  (3.3)  to  replace  the  time  derivatives  with  spatial  derivatives  and  denoting  by 
Uh  the  solution  after  n  time  steps  At,  we  arrive  at  the  final,  fully  discrete  problem  in  the 
form 

+  UX*'=(^I-iAtA,-(l-a}^Al\  Ut  (3.7) 

Remarks; 


1.  Equation  (3.7)  does  not  require  complex  algebra  computations.  The  complex  setting 
is  necessary  to  explore  the  self-adjointness  of  operators  A.  Ah  in  formal  analysis, 
but  does  not  enter  the  computations  (contrary  to  frequency  domain  formulations,  for 
instance). 
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2.  Due  to  the  self-adjointness  of  Ah,  a  stability  analysis  is  immediately  available.  If 
Ai, . . . ,  are  the  (real)  eigenvalues  of  Ah  and  y>i, . . . ,  are  the  corresponding  eigen¬ 
vectors  (note  that  Ah  has  the  same  structure  as  operator  A  on  bounded  domains,  i.e., 
the  eigenvalues  appear  in  pairs  and  the  corresponding  eigenvectors  are  complex  conju¬ 
gates  of  each  other;  see  the  discussion  in  the  previous  section)  then  (3.7)  admits  the 
usual  orthogonal  decomposition  and  for  C/J  =  ipj,  we  have  where 


H  = 


l-f(Aa,)-^^(Aa,)^ 

i  +  f(Aa,)^ 


(3.8) 


Thus,  even  for  an  arbitrary  FE  mesh  (particularly  unstructured  grids)  it  is  possible  to 
express  the  eigenvalues  of  the  transient  operator  in  terms  of  the  eigenvalues  of  operator 
A,. 


3.  A  simple  calculation  reveals  that 


l^il  { 


<  1 

for 

i  <  q;  <  1 

2 

=  1 

for 

1 

Q!  =  X 

2 

>  1 

for 

0<a<i 

(3.9) 


Thus  the  method  is  unconditionally  stable  but  dissipative  for  a  >  ^,  stable  and  energy 
conserving  (like  the  exact  solution)  for  «  =  |  and  unconditionally  unstable  for  o  <  |. 


4.  Note  finally  that 


=  {PhAPhAUh,Vh) 
=  {APhAUh,Vh) 

=  {PhAUh,AVh) 

7^  {AUh,AVh) 


(3.10) 


which  implies  that  the  solution  of  (3.7)  in  just  one  step  would  result  in  a  fully  populated 
stiffness  matrix  (due  to  the  presence  of  the  projection).  A  practical  way  to  avoid  such 
a  situation  is  to  use  the  decomposition 


A/2 

I  +  cy  —  A'i 


I  —  i  \/—  At  Ah 


(3.1 


II 


and  obtain  in  two  steps: 

uT’^  =  (/-;/!  AtA,y'  (l  -  iAtA,  -  (1  -  a)~  a{\  Ul 
Ul*'  =  [l  +  i^AtA,y'ul*i 


(3.12) 


A  different  situation  arises  when  the  method  of  discretization  in  time  is  used.  Using  the 
same  discretization  in  time,  we  arrive  first  at  the  infinite-dimensional  counterpart  of  (3.7), 

^7  +  a  ^  -  (1  -  a)  ^  L/"  (3.13) 

Note  that  the  use  of  (3.13)  requires  increased  regularity  of  the  initial  data  L7°  and  consecutive 
step  solutions  f/",  as  they  must  be  elements  of  D{A^).  This  stringent  assumption  is  avoided 
if  we  replace  (3.13)  with  the  relaxed  variational  version  of  the  form 

(t7"+Sv)+a^  (a17"+\AV) 

=  (C7",  V)  -  iAt  (AU^,  V)  -  (1  -  a) 

V  V  €  D{A) 

as  and  now  must  belong  only  to  D{A). 

If  this  relaxed  formulation  (3.14)  is  adopted  now  as 
replace  simply  U  and  V  with  their  finite  dimensional 
at 

{ur\V,)+a^  {AUl^\AV,) 

=  {Ul,  Vh)  -  iAt  {AUl,  n)  -  (1  -  a) 

VV,  e  X, 

One  easily  recognizes  this  as  the  well  known  Taylor-Galerkin  method  (see,  e.g.,  [3]). 

Obviously,  equations  (3.15)  and  (3.7),  or  equivalently  (3.12),  are  different.  The  key 
point  is  the  relaxation  procedure  performed  at  the  infinite-dimensional  level  in  the  Taylor- 
Galerkin  method  formulation.  Let  us  mention  in  particular  (see  [3]  for  proof)  that  (3.15)  is 
unconditionally  stable  for  a  >  |  and  conditionally  stable  a  <  |.  For  none  of  these  rases 


{AU'^,AV) 


(3.14) 


a  basis  for  spatial  approximation,  we 
counterparts  U hi  V h  S  Xh  arriving 


Ai^ 


{AUI,AVh)  (3-15) 
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is  the  method  energy  conserving.  At  the  infinite-dimensional  level  the  energy  conservation 
takes  place  for  q  =  |.  But  even  for  this  value  of  q,  the  relaxed  formulation  combined  with 
the  spatial  approximation  results  in  extra  “diffusion”  and  eigenvalues  decay  in  modulus  as 
the  “wave  number”  increases. 

Finally,  we  emphasize  that  while  in  the  method  of  lines  the  eigenvalues  of  the  discrete 
transient  operator  were  directly  related  to  the  eigenvalues  of  the  discrete  operator  Ah,  the 
eigenvalues  of  the  transient  operator  resulting  from  (3.15)  can  only  be  related  to  the  spectrum 
of  transient  operator  corresponding  to  (3.14)  which  in  turn  can  be  represented  in  terms  of 
spectrum  of  the  orginal  operator  A. 

This  example  suggests  that  as  long  as  no  relaxation  predure  is  involved  in  the  method 
of  discretization  in  time,  both  procedures  lead  to  the  same  discrete  scheme.  The  class  of 
implicit  Runge-Kutta  methods  we  address  in  the  next  section  satisfies  this  assumption. 

As  a  final  remark,  we  note  that  when  using  the  method  of  lines,  other  approximations  of 
operator  A  are  possible.  In  particular,  upwinded  methods  will  result  in  composing  A  with 
other,  different  projection  operators.  We  do  not  address  this  subject  in  this  paper. 


4  Implicit  Runge-Kutta  Methods.  Rational  Approx¬ 
imations  of  the  Exponential  Function 

Consider  a  general  system  of  ordinary  differential  equations  of  the  form 

y  =  f{y)  (4.1) 

where  y  is  a  1  x  W  vector  of  unknowns  and  /  is,  in  general,  a  nonlinear  function  of  y.  The 
implicit  Runge-Kutta  methods  for  solving  (4.1)  can  be  classified  as  one-step  methods  taking 
an  n-th  step  solution  ?/"  into  through  5  auxiliary  stages  Zj,  j  =  1,2, ...  ,s,  called  the 
internal  approximations,  in  the  following  way 

Z.  =  y" -h  Af^a,j/(Zj),  z  =  l,2,...,s 

j=i 

y"+>  =  y-  +  At^bJ{Z,) 

j=i 

where  a,j,  bj  E  R,i,j  =  1, . . . ,  s.  In  general,  determining  y”'*''  thus  reduces  to  the  solution  of 
N  X  s  nonlinear  (or  linear  if  the  original  system  is  linear)  simultaneous  algebraic  equations. 
There  are  several  good  reasons  for  paying  such  a  price,  and  high  acctiracy  aiul  good  stability 
properties  represent  at  least  two  of  them.  A  complete  discussion  of  the  method  in  the  context 
of  ODEs  can  be  found  in  [1]. 
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Special  cases  are  of  interest.  If  Oij  =  0  for  j  >  t,  system  (4.2)  can  be  decoupled  and  re¬ 
duces  to  the  solution  of  s  consecutive  single  equations  for  internal  approximations  Z\, . . . ,  Z 
Schemes  belonging  to  this  subclass  are  known  as  semi-implicit  methods. 

Another  closely  related  subclass  of  implicit  Runge-Kutta  methods  of  interest  are  the 
so-called  singly-implicit  methods,  where  the  matrix  of  coefficients  aij  is  similar  to  a  lower 
triangular  one.  Reference  to  this  class  of  schemes  as  singly-implicit  refers  to  the  fact  that  the 
coefficient  matrix  has  only  a  single  real  s-fold  eigenvalue.  In  implementation,  methods  of 
this  type  are  as  economical  as  the  semi-implicit  schemes. 


Applying  the  approximation  (4.2)  to  the  Cauchy  problem  of  linear  acoustics,  we  arrive 
at  the  one-step  calculations  in  the  form 


Zfc- =  [/",  A:  =  l,2,...,s 

j=i 

[/n+l  ^  +  At^bjiAZj 

j=i 


(4.3) 


where  Zj,  j  =  l,2,...,s  are  the  intermediate  internal  approximations.  As  the  variational 
formulation  of  (4.3)  does  not  involve  relaxation,  both  the  classical  method  of  lines  and  the 
method  of  discretization  in  time  produce  the  same  result  and  the  spectral  decomposition  of 
the  evolution  operator  is  available  both  at  infinite-  and  finite-dimensional  levels.  Substitut¬ 
ing  for  L7"  an  eigenvector  ^  of  A  and  denoting  by  A  a  corresponding  eigenvalue  of  operator 
A,  wo  arrive  at  the  following  system  of  equations  for  coefficients  dj,j  =  1, . . .  ,s 


dk  -  i(XAi)'^akjdj  =  1  k=l,2,...,s  (4.4) 

j=i 


when-  Zk  —  dk^,  resulting  in  the  formula  for  the  growth  factor  E  in  the  form 


E=l+  i{\At)  bjdj  (4.5) 

}=i 


.Solving  (4.4)  for  dj.j  =  l,...,.s,  we  can  represent  them  as  rational  functions  of  {—iXAt) 
with  ■'ocfficients  depending  upon  the  choice  of  matrix  a*.^. 

Consequently,  the  growth  factor  E  can  be  represented  as  a  rational  function  of  AA/,  too 


E{~7\t) 


N{-iXAt) 

Z^(-^AA/) 


(1.6) 


where  N  and  D  are  j)olynomials.  Note  that  in  the  case  of  (xplicit  Runge-Kutta  metliods. 
(  1.6)  reduces  to  a  polynomial  repre.sentation. 
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Thus,  all  implicit  Runge-Kutta  methods,  when  applied  to  the  acoustics  equations,  fall 
into  the  category  of  one-step  methods  corresponding  to  (4.6),  resulting  in  the  transient 
operator  of  the  form 

=  D(-2AM)-‘  N{-iAtA)U^  (4.7) 

As  polynomial  D  can  always  be  decomposed  into  a  product  of  linear  factors  (possibly  with 
complex  coefficients)  the  solution  of  (4.7)  can  always  be  reduced  to  the  solution  of  .s  de¬ 
coupled  single  equations.  Thus,  for  the  problems  of  interest,  with  the  proper  factorization 
of  polynomial  £),  all  implicit  Runge-Kutta  schemes  can  be  reduced  to  the  solution  of  s 
consecutive  systems. 

Finally,  we  note  that  combining  some  linear  factors  in  D  into  quadratic  terms  and  using 
the  method  of  discretization  in  time,  one  can  arrive  at  various  relaxed  schemes  of  the  Taylor- 
Galerkin  type  discussed  in  the  previous  section. 

We  proceed  now  with  the  discussion  of  some  particular  cases. 

Fade  Approximations 

Figure  2  reproduces  (from  [1])  a  partial  list  of  Fade  polynomial  approximations  to  the  ex¬ 
ponential  function  exp(x).  Each  of  the  approximations  may  correspond  to  a  particular  time 
discretization  scheme  of  type  (4.7).  We  can  classify  them  into  three  groups: 

(a)  Schemes  corresponding  to  the  terms  above  the  main  diagonal  (Gauss-Legcndre 
methods).  They  are  at  most  conditionally  stable. 

(b)  Schemes  corresponding  to  the  main  diagonal.  They  are  all  unconditionally  stable 
and  energy  conserving  as  the  fractions  are  all  in  modulus  equal  one.  In  particular, 
the  second  term  on  the  diagonal  corresponds  to  the  Crank-Nicolson  method  of 
second  order  and  the  third  to  the  two-stage  Gauss-Legendre  method  of  fourth 
order. 

(c)  Schemes  corresponding  to  the  terms  below  the  diagonal.  It  can  be  shown  ([1], 
p.  243)  that  all  terms  from  the  first  two  subdiagonals  are  bounded  in  modulus 
by  one  for  all  x  for  Re  x  <  0.  Schemes  of  this  type  are  called  A-stable.  (Those 
corresponding  to  the  first  subdiagonal  are  known  as  Gauss-Radau  and  to  the 
second  one  as  Gauss-Lobatto  schemes,  respectively.)  In  our  case,  x  is  purely 
imaginary  and,  therefore,  ail  these  schemes  are  unconditionally  stable.  They 
are  only  slightly  dissipative.  Schemes  corresponding  to  the  third  and  further 
subdiagonals  are  not  unconditionally  stable  ([1],  p.  245). 

.\ote  that  for  all  these  schemes  relaxes  versions  arc  possible. 
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Restricted  Fade  Approximations 


The  function 


r{z)  = 


P{zh) 


(4.8) 


(1  -  z/^y 

where  F  is  a  polynomial  of  degree  m,  is  called  a  restricted  Fade  approximation  to  the 
exponential  function  if 

|r(2)-exp(2)i  =  0(2-+>)  (4.9) 


In  practice,  m  =  s  —  1,  s,  5  -f  1.  Schemes  corresponding  to  (4.S)  are  called  singly-implicit. 
They  are  especially  attractive  as  the  solution  of  one  step  requires  only  one  matrix  inversion 
and  s  —  1  resolutions.  It  can  be  proved  ([1],  p.  247)  that 


where  T,  are  the  Laguerre  polynomials  of  order  s.  The  choice  of  pole  7  is  somehow  arbitrary. 
One  can  show,  for  example,  that  for  m  =  s  and  s  =  1,2, . . . ,  6  and  8  the  pole  7  can  be  selected 
in  such  a  way  as  to  render  the  unconditional  stability  (equivalent  to  A  stability  in  this  case); 
for  s  =  7,9,10,...,  such  a  choice  is  impossible.  For  a  summary  of  results,  we  refer  once 
again  to  [1]. 


Numerical  Examples 

VVe  conclude  this  section  with  a  series  of  experiments  comparing  the  behavior  of  various  time 
discretization  schemes  falling  into  the  categories  discussed  on  a  model  problem  depicted  in 
Fig.  3. 

Each  of  the  time-integration  schemes  is  combined  with  a  particular  uniform  h- 
p  approximation  in  space  variables  (see  [2]  for  details).  For  each  of  the  cases,  the  calculated 
pressure  component  is  plotted  at  time  t  =  2.0.  Additional  graphs  provide  an  insight  into 
dissipative  properties  of  a  scheme,  showing  the  variation  of  the  growth  factor  and  phase 
as  functions  of  (the  range  from  0  to  tt,  corresponding  to  the  classical  von  Neumann 
analysis,  is  assumed).  For  comparison,  the  exact  pressure  component  solution  is  shown  in 
Fig.  4. 

The  following  abbreviations  are  used  throughout  the  text 
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Figure  3:  Test  problem  definition  (geometry,  boundary,  and  initial  conditions). 
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SIRK  (5,m) 

GAUSS-LEGENDRE  (5,m) 
RADAU  (5,m) 


LOBATTO  (s,m) 


s-stage,  m-th  order  Singly  Implicit  Runge-Kutta 
Method; 

s-stage,  m-th  order  Implicit  Runge-Kutta  Method 
based  on  Gauss-Legendre  quadrature; 

s-stage,  m-th  order  Implicit  Runge-Kutta  Method 
based  on  Gauss- Radau  quadrature  (designated  as 
Radau  lA  in  [1]); 

s-stage,  m-th  order  Implicit  Runge-Kutta  Method 
based  on  Gauss-Lobatto  quadrature  (designated 
as  Lobatto  IIIC  in  [1]). 


Test  1 — comparison  of  various  schemes  on  a  uniform  mesh  of  quadratic  elements. 


The  first  test  compares  performance  of  various  time  stepping  schemes  with  a  fixed  time 
step  Af  =  0.0625  combined  with  a  fixed  uniform  mesh  of  quadratic  elements  shown  in  Fig. 
5.  (Element  size  h  —  0.125,  order  of  approximation  p  =  2.)  Note  that  the  combination 
At  —  h  corresponds  to  two  time  steps  per  element 

Several  schemes,  including  CRANK-NICHOLSON,  TAYLOR-GALERKIN  and  different 
versions  of  SIRK,  GAUSS-LEGENDRE,  RADAU  and  LOBATTO  were  investigated.  Some 
of  the  results  are  depicted  in  Figs.  6-13. 

The  following  observations  can  be  made: 

1.  Energy  dissipating  schemes  (RADAU,  LOBATTO,  SIRK,  TAYLOR-GALERKIN)  give 
less  oscillatory  results  than  the  energy  conserving  schemes  (CRANK-NICHOLSON. 
GAUSS-LEGENDRE),  and,  therefore,  should  be  preferred. 

2.  The  energy  dissipating  schemes  tested  (RADAU,  LOBATTO,  SIRK)  give  comparable 
results.  Consequently,  SIRK  schemes  should  be  preferred  because  of  their  low  cost 
compared  with  the  RADAU  and  LOBATTO  schemes. 

3.  As  expected,  it  is  inappropriate  to  combine  high  order  time  stepping  schemes  with  low 
order  spatial  approximation.  Increasing  the  order  of  approximation  in  time  does  not 
improve  the  results  and,  e.g.,  SIRK  (5,5)  gives  even  more  oscillatory  results  than  SIRK 

(4,4). 

Finall}",  we  mention  that  several  schemes  belonging  to  the  cla.ss  of  I.inear  Mnltistep  Meth¬ 
ods  have  also  been  tested,  including  the  Method  Based  on  Five-Eight  Rules,  the  .■\dnms- 
Moulton  Method,  and  Gear's  Backward  Difference  Schemes  (see  [1]  for  details).  Of  those,  all 
schemes  of  order  greater  than  2  proved  to  be  unstable  when  applied  to  the  test  prol)lem. 
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Figure  4:  Test  problem.  E.xact  solution  at  t  =  2.0. 
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igurc  5:  Test  1.  An  h-p  FR  mesh,  h  =  0.125.  p  =  2. 
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Figuif'  7:  TAYLOR-CALERKIN  scheme,  a  =  ~,  A/  =  0.0625,  h 


:  Test  1.  IIADAU  (2,3)  scheme,  At  =  0.0G25,  h  =  0.125.  p  =  2. 
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Test  2-Comparison  of  SIRK  schemes. 


As  the  SIRK  schemes  were  judged  to  yield  the  most  promising  results  in  the  first  test, 
the  objective  of  the  second  test  was  to  compare  the  performance  of  SIRK  (.s.  .s  4- 1 )  with  that 
of  SIRK  (s,s).  This  time  a  uniform  FE  mesh  of  cubic  elements,  shown  in  Fig.  M.  was  used. 
The  results  for  SIRK  (3,4)  and  SIRK  (4,4)  are  shown  in  Fig.  15  and  Fig.  16.  respectively. 
As  can  be  seen,  SIRK  (4,4)  produces  better  results  than  SIRK  (3.4). 

Test  3 — Comparison  of  SIRK  (s,s)  schemes  combined  with  different  p-order  FE  meshes. 

In  the  last  test,  an  attempt  was  made  to  find  an  optimal  combination  of  the  order  of 
time  stepping  SIRK  (s,5)  scheme  and  the  spectral  order  p  of  spatial  FE  approximation.  .A 
uniform  FE  mesh  with  h  =  0.125  and  various  orders  of  appro.ximation  p  =  4, 5, 6,  7  was  used. 
The  results  for  various  combinations  of  SIRK  (5,s)  schemes  and  different  p’s  are  shown  in 
Figs.  17  to  25.  It  seems  that  for  the  optimal  combination  one  should  select  p  =  .s  —  1.  .As 
expected,  the  test  results  are  obtained  with  the  highest  order  of  approximation,  namely  for 
SIRK  (8,8)  combined  with  p  =  7. 

5  Adaptivity 

VV’e  conclude  this  paper  with  some  initial  results  on  adaptivity  using  the  Taylor-Galerkin 
scheme  with  cx  —  ^  discussed  in  Section  3. 

Error  Estimation 

The  temporal  discretization  error  has  been  neglected  whereas  the  error  due  to  the  spatial 
approximation  is  estimated  using  the  element  residual  method  in  the  form  (comp.  [8,  9,]) 

||E7/,,p  —  p.).j  11^  ~  (5.1) 

where 


U k.p  is  the  approximate  solution  at  some 


- j - - - ^ 


•  II  •  lU'.n  is  flic  global  “energy’’  norm  defined  as 

\\U\\l,,  =  B{U,U) 
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project:  schock  MESH 


adapt  h-p/2d 


D.O.F-  1S85 


Figure  11: 


Test  2.  An  h-p  FE  mesh,  h 


0.0625,  p 


3. 
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0.  5| 


(5,5) 


DT  *  LAMBDA 


I 


3  .  0 


DT  *  LAMBDA 


where  5  is  a  bilinear  form  defining  the  left-hand  side  of  a  variational  formulation 
corresponding  to  a  particular  time-discretization  scheme. 


•  rjx  is  the  element  error  indicator  function  evaluated  as 

=  Bk{<Pk,<Pk)  (5-3) 

with  the  element  contribution  to  the  global  bilinear  form  and  the  the  clement 
indicator  function  which  is  the  solution  to  the  local  problem 

Find  (pf(  €  X°p^i{K)  such  that 

^  ’  (5.4) 

I  W)  =  R{W)  V  W  €  A'“„,(/C) 

where  R  is  an  appropriate  residual  corresponding  to  element  K  and  A’°p^j(A')  the 
space  of  “bubble”  functions. 


For  all  details  we  refer  to  []. 


While  the  mesh  refinements  are  done  on  the  basis  of  the  element  residual  method,  the 
unrefinements  are  based  on  elements  contribution  to  the  global,  “physical”  ene'^gy,  i.e.,  simply 
the  L^-norm 


il-  =  I  V^u  dx 


(5.5) 


Adaptive  Strategy 

In  the  following  r/  and  ^  denote  the  maximum  values  of  elements  error  and  energy  indicators 

T]  =  maxTjK  ,  ^  =  max^A'  (5.6) 

A  A 

and  rjad,  ^ad  some  fixed  parameters.  The  following  adaptive  strategy  has  been  used  when 
solving  for  a  typical  next  step  solution 

1.  Solve  for 

mesh  refinements 

2.  Estimate  error 

3.  Check  if  t/a-  <  rjad  for  every  clement  K 
if  “yes”  then 

go  to  step  4 
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else 

refine  all  elements  K  for  which  rj^  >  rjad 
go  to  step  1 

endif 


mesh  unrefinements 

4.  Unrefine  all  elements  for  which  <  ^ad 

5.  Set  U”'  =  I/"'*'*  and  go  to  step  1 

Numerical  Experiments — The  Vibrating  Cylinder  Problem 
The  problem  is  defined  as  follows  (see  Fig.  26). 

1.  Governing  equations  (2.6)  are  to  be  solved  in  the  domain 

fl  =  \  {(r,0)  :  r  <  a}  (5.7) 

2.  Boundary  condition 

Ur,  =  A  sin  ivt[Hit)  -  H{t  -  T/2)]  (5.8) 

where  A,  cu,  T  are  constants  and  H{-)  is  the  Heaviside  function. 

3.  Sommerfeld  radiation  condition  at  r  =  oo 

4.  Initial  condition 

Uo  =  0  (5.9) 

For  computations  we  accept  the  computational  domain 

n  =  {(r,  0)  :  a  <  r  <  Too  ,  —a  <  0  <  a} 

The  problem  was  solved  using  the  Taylor-Galerkin  method  with  a 
time  step  At  —  0.015625.  The  following  boundary  conditions  were 
calculations:  Boundary  condition  (5.8)  at  r  =  a,  —a  <  0  <  a  and 

u„  =  0  at  a  <  r  <  Too  ,  0  =  ia 

=  0  at  r  =  r^  ,  —a  <  0  <  a 
an 

The  adaptive  strategy  was  based  on  /i- refinements.  Figure  27  shows  an  initial  FE  mesh 
of  quadratic  elements  and  the  consecutive  figures  28,  29,  and  30  present  the  evolution  of 
FE  mesh  at  time  stages  t  =  0.5,  1.0  and  t  =  1.5,  respectively.  Finally,  Fig.  31  shows  the 
computed  pressure  distribution  at  time  t  =  1.5. 


(5.10) 

=  0.5  and  a  constant 
used  in  one  time-step 
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Figure  30:  The  Vibrating  Cylinder  Problem,  /i-refined/unrefined  mesh  at  time  =  1.5. 


Figure  31:  The  Vibrating  Cylinder  Problem.  Pressure  distribution  at  time  =  1.5. 
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6  Conclusions.  Further  Research 


The  paper  presents  some  preliminary  results  toward  designing  fully  automated  self-adaptive 
methods  for  solving  transient  problems  of  type  (2.i0)  with  applications  to  linear  acoustics. 
Two  concepts  of  discretization  have  been  presented  and  discussed:  the  method  of  lines  and 
the  method  of  discretization  of  time.  While  the  latter  one  provides  a  natural  basis  for 
every  step  mesh  refinements/unrefinements,  the  first  allows  for  an  explicit  representation  of 
spectrum  of  the  transient  operator  in  terms  of  the  eigenvalues  of  the  approximation  operator 
Ah,  including  the  case  of  arbitrary,  unstructured  meshes.  Consequently,  all  diffusive  and 
dispersive  properties  of  the  transient  operator  are  known  (in  particular  the  stability  analysis 
is  easily  available)  and  the  whole  analysis  reduces  to  the  investigation  of  the  spectrum  of 
operator  Ah  as  a  function  of  adapting  meshes. 

The  implicit  Runge-Kutta  schemes  are  shown  to  belong  to  the  class  of  methods  for  which 
both  discretization  concepts  yield  the  same  result  and  seem  to  be  a  natural  candidate  for 
constructing  approximations  to  the  considered  problem. 

Several  questions  remain  unresolved  and  are  under  current  investigation.  To  mention  a 
few: 

•  Dependence  of  the  spectrum  of  the  approximate  operator  Ah  on  spatial  approximation, 
including  h-p  meshes.  For  regular,  uniform  grids,  eigenvalues  of  Ah  are  uniformly 
distributed  and  often  are  available  in  a  closed  form  (that  is  why  the  von  Neumann 
stability  analysis  is  available).  Their  structure  for  non-uniform  meshes  is  not  known. 
Does  the  adaptivity,  in  particular  that  causing  local  transitions  in  element  size  {h)  and 
order  (p),  affect  strongly  the  shape  of  the  corresponding  eigenvalues? 

•  If  the  mesh  is  refined /unrefined,  the  spectral  representation  of  Ah  changes  and  the 
corresponding  solution  has  to  be  represented  using  different  eigenfunctions  (transfer 
of  energy  between  different  modes).  Can  the  spectral  properties  of  Ah  be  controlled? 
control  it?  What  effect  does  it  have  on  the  error? 

•  A  combined  a  posteriori  error  estimate  including  both  time  and  space  discretization. 
Some  of  the  SIRK  schemes  allow  for  an  error  control,  provided  one  extra  iteration  is 
performed  (comp.  [1]).  It  may  be  possible  to  combine  those  techniques  with  residual- 
type  estimates  with  respect  to  the  space  variables. 

•  Is  it  possible,  by  using  adaptive  methods,  to  obtain  exponential  convergence  for  the 
transient  hyperbolic  equations? 

We  liopc  to  answer  .some  of  these  cjnestions  in  a  forthcoming  work. 
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